Philippe Poiré

Their great competences and knowledge on the inertial and satellite navigation domain have been essential

Introduction

Background

The upcoming of services based on the location of the user/vehicle are at the origin of the increasing demand for high performance personal navigation. First driven by the issue of directives to provide a location mean to emergency calls (i.e. the E911 mandate set by the Federal Communications Commission in the US and the E112 directive set by the European Commission in Europe), different locations based services (LBS) are being developed for commercial applications. Some of the proposed positioning products include the availability of connecting users to nearby points of interest (such as retail businesses, public facilities, or travel destinations), advising them of current conditions (e.g. traffic), or providing routing and tracking services. Among the different growing markets, incar navigation is becoming a leading application either for personal or professional use. Services associated to vehicle navigation such as fleet management systems, road tolling or pay-as-you-drive applications are requiring high performance positioning both in terms of accuracy and reliability.

Global Navigation Satellite Systems (GNSS) are designed to provide position, velocity and timing capabilities to user through out the world. Currently, the only fully operational GNSS is the Global Positioning System (GPS). Over the years, increasingly falling cost of GPS receivers, as well as their size and consumption, has rendered the system attractive for the design of land vehicle navigation systems. However, satellite navigation presents several impairments in degraded signal environments such as urban canyons. Land vehicles typically move near high buildings and dense foliage zones where GPS navigation capabilities are usually jeopardized. Signal interruption is one of the primary reasons which affects the continuity and reliability of the GPS navigation solution. Either the received signals 1 are too weak to be tracked using conventional techniques, either they are totally obstructed (i.e. in non-line-of-sight). Multipath and interferences will entail a severe degradation in the navigation performance.

Self-contained augmentations have the advantage of being independent from the navigation environment. Dead reckoning (DR) navigation systems are a typical example of self-contained systems.

They are not affected by external radio frequency signals and they ensure continuity over the position estimation. Thanks to the efforts done by the industry to miniaturise DR sensors and render them affordable for commercial applications, they are now gaining importance in the personal navigation field. However, they only provide information on the relative movement of the vehicle and they suffer from time-dependent error growth which causes a drift in the solution. Hybrid techniques fusing GPS/DR-system are thus implemented to exploit capabilities from both navigation systems. Their combination not only offers the accuracy and continuity in the solution, but also enhances the reliability of the system. GPS can restrict the DR navigation system error growth over time, and allows for online estimation of the sensor errors. DR sensors can bridge the position estimates when there is no GPS signal reception and limit the impact of severe GPS measurement errors.

In an attempt to overcome limitations dealing with low power received signal high sensitivity GPS receiver technology is used to acquire and track weak GPS signals. Receivers use long coherent integration times to reduce the effect of noise and increase the probability of detecting a specific satellite signal. However, increased measurement noise due to lower signal strength (potentially leading to undesired cross-correlation effects and thus biases) and the high levels of signal reflection (i.e. multipath) prevent the receiver from achieving high performances. This situation prompts the need for a reliability analysis to detect the presence of a defective signal and mitigate its influence in the final position estimation.

Thesis framework

The development of this thesis is aimed at coping with urban canyon phenomena and the problems they entail in the precision and reliability of the navigation solution. It is important to understand that due to the statistical properties of the received satellite signals every given solution has an associated error probability distribution function (pdf). Therefore, the main challenge of this study is not only providing accurate position estimations under degraded signal reception scenarios, but being capable of correctly modeling their associated error pdf. In this way, a measurement of trust (i.e. an error bound) can be associated to the navigation solution. The need for reliable bounding reveals critical for several land applications such as the transport of dangerous materials or road tolling. Indeed, for open-sky scenarios, the GPS measurements have proved to follow a zero mean Gaussian distribution which can be easily characterized. However, for urban navigation, the correct determination of the error model describing the received signals is a great challenge and thus the subject of many researches. This thesis will be motivated by two parallel but complementary approaches. Indeed, both aim at enhancing the navigation performance but each of them in a different way: Different approaches exploiting GPS data and on-board vehicle information will be of interest for the development of a robust navigation system suitable for urban environments. It is important

Introduction

to note that all by itself the hybrid GPS/DR approach does not guarantee the consistency of the navigation solution. It surely helps to enhance positioning accuracy (by providing some compensation for GPS erroneous measurements), but it does not directly address the reliability problem. An erroneous GPS measurement with an associated inconsistent error model will still have some impact in the final position solution.

Detection and mitigation of GPS errors introduced by urban phenomena

Signals received in urban canyon scenario will be affected by different errors sources such as multipath or cross-correlation errors. These types of errors present a major problem because their evolution or error model is usually unknown (different from the ionospheric errors, for example, were models are used to calculate the introduced error). A reliability test must be implemented to account for any inconsistency between the received measurement and its error pdf. Usual approaches working on the quality of GPS measurements are based on an exclusion criterion. In other words, if an outlier or erroneous measurement is detected, it is automatically excluded from the navigation filter solution. However, this may present a high cost for urban navigation where GPS satellite visibility is scarce. This thesis is aimed at studying different techniques not to exclude the outlier satellite but to "compensate" its error.

Author contribution

This thesis presents a general study of land vehicle navigation in urban environments. Considering the thesis framework discussed in the previous section the contributions of this thesis can be summarized as follows:

INS and differential odometry coupling approach A high level integration is proposed for three self-contained navigation systems: GPS, INS and wheel speed sensors based differential odometry system. While the coupling between the first two system is a quite standard, the contribution of this thesis lies in the incorporation of the third system. Combination between inertial and odometric measurements is intelligently done to achieve a 2D navigation.

Non linearity impact on the navigation solution The Unscented Kalman Filter (UKF) is tested to provide information on the correctness of the linear assumption made when implementing the standard Extended Kalman Filter (EKF). In this way, conclusions on the impact of nonlinearities for land navigation can be drawn.

Smooth correction criterion

An innovative approach is proposed to exploit the robustness of Doppler measurements to urban canyon interference. A two threshold detection technique is proposed where either the pseudorange, the Doppler or a smoothed measurement (obtained by a special merging between pseudoranges and Doppler measurements) are used. In this way, any lack of consistency between the received pseudorange and its associated error model should be mitigated.

Outliers' error model The contribution relies on the three hypothesis model considered for the characterization of the received signal. The three hypothesis include the nominal case and two other cases accounting for either a mean or variance jump in the noise model. This mainly approach allows the identification of multipath presence either in the presence of line-of-sight or non-line-of-sight direct signal.

Hierarchical multipath mitigation technique

Based on a three hypothesis model (described in the previous point) a promising detection, identification and correction technique is deployed. Indeed, not only the outlier presence is detected, but its source of error is identified, evaluated and further compensated. As no exclusion is done, this approach should be particularly well suited for urban vehicle navigation.

Thesis outline

The thesis is organized in 5 chapters.

Chapter 1 gives a general overview of the GPS principles, structure and functioning. Special attention is paid to the analysis of the different urban phenomena degrading the received signal.

Indeed, this point is of major importance because it motivates all further studies presented in the following chapters. A brief introduction to the high sensitivity GPS receiver and Egnos contributions is provided.

Chapter 2 introduces an alternative navigation systems based on dead reckoning sensors. A first introduction to inertial sensors (i.e. accelerometers and gyros) and wheel speed sensors (WSS) is provided along with their corresponding error model. The Inertial navigation System (INS) mechanization is derived for a 2D navigation. Expressions for the WSS-based navigation system obtained from the application of differential odometry principles are presented.

Chapter 3 deals with the development and analysis of different types of hybrid strategies. In particular, the UKF will be presented as an alternative to the traditional EKF used for navigation purposes. DR solutions presented in the previous chapter will augment the GPS-based navigation filter. A high level GPS/INS/WSS approach is presented and its performances contrasted with a nimonal GPS system. Observations are done on the relevance non-linearities present for vehicle navigation.

Note: concepts, results and conclusions from previous chapters will help to the definition of algorithms presented in the following two chapters.

Chapter 4 is based on Doppler measurements robustness to urban canyon phenomena. They will therefore be proposed as an alternative/complement to traditionally used pseudorange measurements. A smooth correction criterion provides an innovative framework to test the reliability of the received pseudoranges. A merging between pseudoranges and Doppler measurements is proposed as an interesting solution to situation where the presence of an outlier must be ambiguous. Principles for the computation of an adapted positioning bound are given so that performances can be evaluated both in terms of accuracy and model consistency.

Chapter 5 presents two-step detection procedure that aims at classifying outliers according to their associated source of error. Two different situations will be considered in the presence of multipath. These situations correspond to the presence or absence of line of sight signal for the different GPS satellites. Therefore, two kinds of errors are potentially "corrupting" the pseudoranges, modeled as variance changes or mean value jumps in noise measurements. An original multiple model approach is proposed to detect, identify and correct these errors and provide a final consistent solution.

Chapter 1

GPS Satellite Navigation

Introduction

Navigation has been used since ancient times by exploiting the angular measurements from natural stars. However, this strategy depended on sky visibility (i.e. on whether conditions) and provided low accuracy. With the upcoming of radio communication the possibility of developing new navigation system was foreseen. By knowing the position of the station emitting the signal, the user could compute its position. First approaches were conceived for ground based station. However, a trade-off between the accuracy and the coverage of the system was to be considered. High frequency signals provide accurate measurement but need line-of-sight (LOS) users while low frequency signals provide lower precision but enable a wider coverage. These considerations highly jeopardize the feasibility of a global navigation system. The real breakthrough to what would later lead to a sophisticated satellite navigation system was made in the 1960s, when the US navy presented its satellite navigation system known as transit. High frequency signals were transmitted and a wide coverage could be assured.

The system measured the signals' Doppler shift and provided a two-dimension positioning. The time needed to obtain a position fix could be up to 110 minutes depending on the user's latitude.

Therefore, though this system could be used for ship navigation because of its low velocities and known height, it was not suitable for air navigation. To overcome these type of limitations the U.S Global positioning system (GPS) and its analog Russian Global Navigation Satellite System (GLONASS), were developed. In particular, all the studies concerned within this PhD are based on GPS signals.

This chapter is aimed at providing an overview of the Global positioning system structure, functioning, and limitations. Urban navigation challenges will be specially studied. Descriptions and discussions will be mainly based on reference textbooks such as the ones published by Parkinson and Spilker [START_REF] Parkinson | Global Positioning System: theory and applications[END_REF], Kaplan [START_REF] Kaplan | Understanding GPS principles and applications[END_REF] or Misra and Enge [START_REF] Misra | Global Positioning System: Signals, Measurements and Performance[END_REF]. The GPS modernization program (i.e. GPS III program) and its associated modernized signals such as L2C or L5 have not been considered in this thesis. The interested reader is invited to consult [START_REF] Kaplan | Understanding GPS principles and applications[END_REF] for more details.

Satellite Navigation

Global Positioning System

The Navstar (Navigation System by Timing and Ranging) GPS, most commonly referred as simply GPS, is a LOS weather independent, world-wide continuously available Radio Frequency (RF) positioning system, that provides three-dimensional position, velocity and time synchronization capabilities to end-users with an appropriate receiver. The GPS parameters are framed geometrically in the ECEF (earth centered earth fixed) WGS-84 (World Geodetic) world-wide common grid reference system. The position calculation is done by triangulation of distance measurements between the user and the satellite. GPS applies the concept of one way time of arrival (TOA), utilizing satellite transmissions referenced to highly stable atomic clocks on board the satellite and "synchronized" with the reference GPS system time. Satellites transmit continuous waveforms that can be easily related with a time scale. In fact, the received signal is compared with a locally generated copy.

Assuming the receiver is ideally synchronized with the satellite clock, the time difference between the received and generated signals will be the propagation delay. If the transmitted time and the satellite position are known, the distance between the receiver and the satellite (i.e. the range measurement) can be obtained. The satellite position is included in the transmitted ephemeris as well as other parameters to correct the drift of the satellite clock. In this way a 3D position solution is obtained. The user velocity is computed from the Doppler shift caused by the relative motion between the satellite and the receiver. However, in real life, receiver clocks are not stable and can drift significantly from the GPS time. This drift will strongly impact the correct computation of the propagation time. Hence, the GPS measurement will not exactly correspond to the receiver/satellite distance and it will be actually called pseudorange measurement. This fourth unknown parameter will be added to the other three basic ones (corresponding to the three-dimensional receiver position).

Four satellites are thus needed to correctly determine user latitude, longitude, height and receiver clock offset from internal system time. GPS is composed of three different segments, each of them aimed at performing very specific tasks. These segments are related to satellite constellation (Space segment), ground-control/monitoring network (Control segment) and user receiving equipment (User segment).

Space segment

The space segment was originally specified for 24 satellites. The constellation is designed to operate in six Earth-centered orbital planes, 60 • degrees apart, nominally inclined at 55 • to the equator.

Each orbital plane contains four to five satellites orbiting at an altitude of 20183km from the mean surface of the Earth, with a period of one-half of a sidereal day (approx. 11 hours and 58 seconds). However, the system currently employs more satellites than specified in the nominal constellation, and at the time of writing, the GPS constellation consisted of 31 Block II/IIA/IIR satellites (each block corresponds to a different phase of development). With the increased number of satellites, the constellation was changed to a nonuniform arrangement. Thus, the current optimized constellation has up to 7 orbital slots unevenly spaced around each plane. Such an arrangement was shown to improve reliability and availability of the system relative to a uniform system. In this way about 9 satellites are visible from any place at any time, considering an elevation mask of 5 • . Note that under this elevation the degradations suffered by the GPS signal due to its passage through the atmosphere prevent a precise positioning to be guaranteed. The system currently includes two different carrier frequencies to modulate the GPS signal: L1 (at 1575.42 MHz) and L2 (at 1227.60 MHz). Both signals are modulated by spread spectrum codes with an unique pseudo-random (PRN) sequence associated to each satellite. These two emitted GPS signals are thus expressed as follows:

s L1 i (t) = √ 2Ad i (t)C i (t) cos(2πf 1 t + θ 1 ) + √ Ad i (t)P i (t) sin(2πf 1 t + θ 1 ), (1.1)
s L2 i (t) = A 2 d i (t)P i (t) cos(2πf 2 t + θ 2 ), (1.2)
where,

• i = 1, ..., n y , n y denoting the number of visible satellites,

• A is the mean power of the ith emitted signal.

• d i (t) represents the navigation data corresponding to the ith satellite. It mainly provides the means for the user to determine the location of the satellite at the time of signal emission. This message includes information about the satellite position (i.e. its ephemeris), clock behavior with respect to GPS time, system status and ionospheric error models.

• C i (t) denotes the coarse/acquisition (C/A) pseudo random noise (PRN) code for the ith satellite,

• P i (t) denotes the Precision (encrypted) (P(Y)) PRN code for the ith satellite,

• f 1 = 1575.42M Hz and f 2 = 1227.60M Hz denote the carrier frequencies associated to the L1 and L2 signals respectively.

The s L1 signal is modulated using a quaternary phase shift keying (QPSK) digital modulation. The in-phase channel is modulated by the C/A code, while the quadrature channel is modulated by the P(Y) code. The s L2 contains just the P(Y) code with a binary phase shift keying (BPSK) modulation.

The C/A nd P(Y) codes differ in their periods and chipping rate. The higher the chipping rate and the longer the period of the code, the more precise the positioning solution but the slower the time of acquisition. noise. In this way, the signal becomes highly robust to jamming or interferences, with a low probability of being detected. Moreover, thanks to the code properties, each signal can be uniquely identified through a correlation operation. These code properties will be exploited by the receiver to estimate the signal time of arrival. For this purpose, a spreading code local replica of the desired satellite will be generated by the user's receiver and it will be correlated with the received signal. The code properties must be as close as possible to those of an infinite random noise sequence. In other words, it must have a delta-like autocorrelation function centered in 0 and very weak cross-correlation properties to ensure precise synchronization and limited interference from other satellites. However, the code must have at the same time a finite length to enable the receiver to reproduce it and perform the signal acquisition within a short period time. The Gold code family presents an interesting equilibrium for this trade-off and they are consequently used to generate the C/A code.

Control segment

The control segment ensures the surveillance of the GPS constellation. Its principal functions consist of

• keeping the satellite in its adequate orbit and position,

• monitoring the satellite health,

• estimating and predicting the satellite clock drift parameters and ephemeris

•
r(t) =    ny i=1 M i j=0 A j i d i (t -τ j i )c i (t -τ j i ) cos(2π(f 1 + f j d,i )t -θ j i )    + w(t), (1.3)
where

• the index i denotes the ith LOS signals,

• M i is the number of multipaths associated to the ith visible satellite,

• A j i is the signal amplitude corresponding to the jth path of the ith satellite,

• τ j i is the signal propagation delay associated to jth path of the ith satellite,

• f 1 is the L1 frequency,

• f j d,i
is the Doppler frequency associated to jth path of the ith satellite,

• θ j i is the Doppler phase shift corresponding to the jth path of the ith satellite,

• w(t) is an additive white noise.

In the down conversion stage, the received L1 signal is amplified, filtered, downconverted (to an intermediate lower frequency) and digitalized. The signal is then going to be processed by each of the N channels to identify the N possible visible satellites. Each channel will decompose the signal associated to the tracked satellite in its in-phase (I) and quadrature (Q) components. The resultant I-Q signals will be then correlated with a local code replica for the desired satellite. These correlation samples are later going to be used by the code and carrier loops at the acquisition and tracking stage.

In the absence of multipath, the correlator outputs for the I and Q components corresponding to the ith satellite are written as,

I i (τ ) = A i 2 d i sin(π∆f i T b ) π∆f i T b R(τ -τi ) cos(π∆f i T b + (θ i -θi )) + n I,i (t), (1.4)
Q i (τ ) = A i 2 d i sin(π∆f i T b ) π∆f i T b R(τ -τi ) sin(π∆f i T b + (θ i -θi )) + n Q,i (t), (1.5)
where

• d i is the navigation data considering that no data bit transition occurs during the integration time,

• ∆f i is the error induced by the difference between the real and estimated Doppler frequency,

• T b is the integration time,

• R(•) is the C/A code correlation function,

• τi is the estimated propagation delay,

• θi is the estimated Doppler phase shift,

• (n I,i (t), n Q,i (

Types of measurements

Three different GPS measurements can be computed from the received signal: the pseudorange (obtained from the code phase), carrier phase and Doppler measurement. Further details on the characteristics associated to each of these measurements, as well as their domain of application, will be given in chapter 4. In the following, measurement models, and their relation to the vehicle state, are presented.

Pseudorange measurement As already introduced in the previous section, once the received signal is identified, the DLL estimates the shift between the local generated C/A code and received signal code. When the point of maximum correlation is found, synchronization is achieved. The estimated code shift allows to obtain the propagation time between the satellite and the receiver.

Multiplying this time shift by the speed of light results in the apparent satellite-to-user range (i.e.

the pseudorange). The pseudorange model for nominal open-sky scenarios can be written as

ρ i = r i + c ∆t R -∆t S i + eph i + iono i + tropo i + rel i + n ρ i , (1.6)
where

• ρ i is the pseudorange measurement in meters associated to the ith signal,

• r i represents the effective geometric range between the ith satellite and the receiver,

• c is the speed of light,

• ∆t R is the receiver clock error with respect to the GPS time,

• ∆t S i is the ith satellite clock error with respect to the GPS time,

• eph i
is the error related to ephemeris inaccuracies.

• iono i
is the error induced by the ionospheric delay on the ith satellite signal,

• tropo i is the error induced by the tropospheric delay on the ith satellite signal,

• rel i is the relativistic time error for satellite i,

• n ρ i is the noise affecting the ith pseudorange measurement.

The compensation of errors introduced by the ionosphere, troposphere, satellite clock and inaccuracies in the ephemeris parameters is addressed in section 1.2.2. These errors are going to be addressed in the following as "modeled" errors because they can be locally modeled and partially compensated.

The induced relativistic effect is compensated thanks to the correction model parameters included in the navigation message. Hence, under clear-sky conditions, the pseudorange model to be considered hereafter for the receiver position computation is expressed as

ρ i = (x -x s i ) 2 + (y -y s i ) 2 + (z -z s i ) 2 + b t + n ρ i , (1.7)
where

• ρ i is the pseudorange measurement after the modeled errors have been compensated,

• (x -x s i ) 2 + (y -y s i ) 2 + (z -z s i )
2 denotes the range measurement r i where both the satellite and receiver position are made explicit,

• (x, y, z) denotes the receiver coordinates in the geodetic Earth-centered-Earth-fixed (ECEF) frame according to the World Geodetic System WGS-84,

• (x s i , y s i , z s i ) denotes the satellite coordinates in the geodetic ECEF for the ith GPS satellite,

• b t is the receiver clock bias, where b t = c∆t R ,

• n ρ i represents the noise affecting the measurement plus the residual errors after correction.

Therefore a minimum of 4 satellites is needed to solve for the 3 dimension vehicle position and receiver clock bias. For more details on the different existing coordinate frames please refer to chapter 2.

Doppler measurement Doppler measurements are based on an existing relative motion between two sources. For a moving emitter (namely a satellite) or a moving receiver, the received frequency is shifted with respect to the emitted one. By measuring this frequency shift and knowing the satellite velocity (using ephemeris data) and clock drift, the receiver's velocity can be determined. Considering a similar model to (1.7), where the "modeled" drift errors have already been compensated, the Doppler measurement can be written as

D i = - L1 c v e -v S i T • l LOS i + d t + n D i (1.8)
where

• D i is the Doppler measurement in Hz for the ith satellite

• v e is the receiver velocity vector in geodetic ECEF coordinates,

• v S i is the satellite velocity vector in geodetic ECEF coordinates,

• the operator • denotes the inner product,

• l LOS i is the unit line of sight vector between the i th satellite and the receiver,

• d t is the receiver clock drift error with respect to the GPS time so d t = ∆ ṫR ,

• n D i is the Doppler measurement noise. . Indeed, the Doppler shift is the projection of the relative satellite-to-user motion into the LOS vector. This unit LOS vector is obtained as:

x -x s i r i , y -y s i r i , z -z s i r i T , (1.9)
where

r i = (x -x s i ) 2 + (y -y s i ) 2 + (z -z s i ) 2 .
Note that a previous position estimation is necessary to calculate the l LOS i vector. A minimum of four Doppler measurements is needed to estimate the 3 dimension vehicle velocities and the receiver clock drift. It is interesting to consider that as the vehicle coordinates are indirectly used in the Doppler measurement model, positioning could be done by using the Doppler measurements. However, this approach is substantially less precise than the pseudorange positioning [START_REF] Axelrad | GPS navigation algorithms[END_REF].

Carrier measurement

The carrier phase is an indirect and ambiguous measurement of the userto-satellite range. If the carrier signal is constantly tracked, the accumulation of the time varying carrier phases leads to a knowledge on the changes in the observed range. This measurement is commonly used in high-accuracy applications such as static or dynamic surveying or attitude determination. Considering a previous compensation of the atmospheric effects and satellite clock bias, the measurement is expressed as follows

φ i = (x -x s i ) 2 + (y -y s i ) 2 + (z -z s i ) 2 + N i λ + b t + n φ i , (1.10)
where

• φ i is the accumulated delta range measurement corresponding to the ith satellite,

• N i is the carrier phase ambiguity corresponding to the signal from the ith satellite,

• λ is the L1 carrier phase wavelength,

• n φ i is the noise affecting the delta range measurement. This measurement will not be used within this PhD. Further discussions are provided in chapter 4.

Measurement errors in open-sky scenario

Several different errors affect the received signal as already introduced in (1.6). They can be classified in three categories: satellite-based errors, signal propagation errors and receiver-based errors. The satellite-based errors include satellite clock and ephemerids errors. Among the signal propagation errors, ionospheric and tropospheric errors are found. The receiver errors are mainly associated to receiver noise affecting the correct functioning of the tracking loops. The ability to obtain accurate and reliable tracking capabilities, resides in the predictability, controllability, and detectability of the measurement errors. In the following, the different sources of errors affecting measurements in an open-sky scenario are described.

Satellite clock error

Although the atomic clocks contained by the GPS satellites are highly stable, they are also affected by a small drift. These errors are due to offsets in the clock frequency of each satellite with respect to the reference clock. The satellite errors are modeled via a polynomial, whose coefficients are transmitted as a part of the navigation message. The satellite clock error is usually less than 1 ms and, after implementing the broadcast correction, the remaining error is in the order of 8 to 10 ns (2 to 3 m). The relativistic correction is also included in the polynomial clock correction model to account for the fact that the slight eccentricity of the satellite orbit causes it to travel through different levels of gravitational potential. For instance, when the satellite is at perigee (i.e. closest to the Earth) the gravitational potential and the satellite velocity are higher, causing the clock to run slower. The contrary effect is observed when the satellite is at apogee.

Ephemeris error

The ephemeris error (also called orbital errors) results when the transmitted ephemerids do not accurately describe the true satellite position. Therefore, the accuracy of the model for predicting the satellite position is limited. A typical ephemeris error is of approximately 1 meter on the range measurement.

Ionospheric error

The ionosphere is a dispersive medium which extends from about 50 to 1000 km above the Earth and is characterized by an abundance of free electrons and ions. When the signal travels through the ionized air, the signal speed decreases from the vacuum speed of light and therefore results in phase advance and code delay. Ionospheric effects show diurnal and seasonal variations.

They depend on the solar cycle and the geographic location of the receiver. The ionospheric delay is also a function of frequency. Dual frequency GPS users can utilize this property to correct the error in range and Doppler measurements. However this approach is out of the scope of this thesis.

Single frequency users can compensate for ionospheric error using the approximate model (e.g. the Klobuchar model [START_REF] Klobuchar | Ionospheric effect on GPS[END_REF]), which is based on broadcasted parameters included in the GPS navigation message. This algorithm removes nearly the 50% of the ionospheric error. The ionospheric error is typically around 2 to 20 m [START_REF] Misra | Global Positioning System: Signals, Measurements and Performance[END_REF] depending on the satellite elevation angle with respect to the receiver.

Tropospheric error

The troposphere is the lower of the Earth atmosphere and affects GPS signals at altitudes up to 50-70 km above the surface. The introduced error will be characterized by the humidity, pressure and temperature along the signal path in the troposphere. It has to be computed by the receiver according to average metereological parameters. Tropospheric delay is subdivided into dry and wet components due to distinct influences on RF signal propagation. The dry portion typically contributes 80-90% to the entire effect and it can be easily predicted. By contrast, the wet term constitutes the remaining 10-20% of the total error, and it is far more difficult to estimate because vapour density varies with the local weather. One simple model that is often used for the tropospheric error compensation is the Hopfield model [START_REF] Hopfield | Two quadratic tropospheric refractivity profile for correction satellite data[END_REF]. For most users and circumstances, a simple model should be effectively accurate to about 1m [START_REF] Spilker | Tropospheric effect on GPS[END_REF].

Tracking Loop Jitter The presence of noise in the tracking process will affect the phase shift and code delay measurements, and therefore the final pseudorange or Doppler measurement. For these measurements noise increases as signal power decreases due to a growing thermal noise jitter and dynamic stress error in the carrier and code tracking loop. Pseudorange measurement noise depends on correlators' spacing, precorrelation bandwidth, and on the bandwidth of the delay lock loop (DLL) used in code tracking. It will also be determined by the type of code discriminator used [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF].

The expression for the noise variance σ 2 ρ i affecting the code measurement considering a DLL with dot-product discriminator can be written as follows [START_REF] Van Dierendonck | GPS receivers[END_REF]:

σ 2 ρ i ≈ B l C s 2C/N 0 1 + 1 C/N 0 T c (1.11) where • σ 2 ρ i is expressed in (chips) 2 ,
• B l is the loop noise bandwidth,

• C/N 0 carrier to noise power [Hz],

• C s is the spacing between correlators,

• T c is the prediction integration time.

Similarly, the Doppler measurement noise mainly depends on the thermal noise of the frequency tracking loop (considering that receiver effectively obtains this measurement from the FLL) and on the correlation bandwidth. The expression for the noise variance σ 2 D i affecting the Doppler measurement considering a FLL with cross-product discriminator can be written as follows [START_REF] Natali | AFC tracking algorithms[END_REF]:

σ 2 D i ≈ 1 π 2 B l B 2 I C/N 0 1 + B I 2C/N 0 . (1.12)
where σ 2 D i is expressed in Hz and B I is the predetection noise bandwidth. The receiver noise is generally 0.5m for code measurements and around 0.5Hz for frequency measurement. Table 1.2 presents the pseudorange error characterization in terms of their 1σ standard deviation according to [START_REF] Parkinson | Global Positioning System: theory and applications[END_REF]. Since the error sources are reasonably independent, the final measurement noise is considered to be modeled as a zero mean Gaussian distribution with an associated variance given by σ U ERE (where U ERE stands for User Equivalent Range Error).

Class

Error • Fast corrections: they are used to correct the fast changing errors which were usually caused by the satellite clock error due to degradation through selective availability techniques (now no longer existing). Therefore this term was not used in this thesis.

• Slow corrections : they are aimed at correcting the slowly varying errors. Among these errors we can cite the satellite position errors caused by errors in the ephemerids transmission and the degradation of the satellite position calculation with the time. They also complement the satellite clock error model provided by the GPS message.

• Ionospheric corrections: they provide a more accurate ionospheric correction than the one computed from the GPS message. A more complex ionospheric grid model is used.

In the presence of EGNOS corrections, the pseudorange measurement corresponding to the model given in (1.7) will be obtained as,

Y i = ρ i + RC i f ast + RC i iono + RC i tropo + RC i clock , (1.13)
where

• ρ i is the pseudorange measurement computed after applying the satellite position corrections present in the EGNOS message.

• RC f ast are the fast corrections for the ith satellite computed from the EGNOS message,

• RC iono are the ionospheric corrections for the ith satellite computed from the EGNOS message

• RC tropo are the tropospheric correction corresponding to the ith satellite computed by the receiver,

• RC clock are the satellite clock corrections corresponding to the ith satellite computed from the EGNOS message (they belong to the slow corrections).

It is important to notice that the tropospheric correction will be computed at the receiver. The RC tropo term does not depend on the reception of any correction information transmitted in the EGNOS signal. They depend on metereological conditions whose parameters are interpolated from already tabulated tables. The expression given in (1.13) corresponds to the error correction model

given in [START_REF]PEGASUS: Technical notes on SBAS[END_REF]. Furthermore, EGNOS provides not only necessary parameters for error correction, but also an upper bound on the residual error after corrections were applied. This bound is given in the form of a Gaussian variance. Hence, the actual variance of the pseudorange measurement is going to be calculated according to the model given in [START_REF]PEGASUS: Technical notes on SBAS[END_REF],

σ 2 i = σ 2 i,f lt + σ 2 i,iono + σ 2 i,tropo + σ 2 i,rec (1.14)
where

• σ 2 i is the total pseudorange measurement variance

• σ 2 i,f lt is the variance of the residual error after the fast and slow corrections were applied,

• σ 2 i,iono is the variance of the residual error after the ionospheric corrections were applied,

• σ 2 i,tropo is the variance of the residual error after the tropospheric corrections were applied,

• σ 2 i,rec is the variance for the receiver errors (it depends on the receiver characteristics (1.11)). The erroneous tracking of a multipath is also possible when either the received signal is too weak or inexistent (i.e. an obstacle was totally obstructing the LOS path). This section is aimed at providing some details on these types of error sources. Indeed, it must be pointed out that due to their strong environmental dependency, urban navigation errors cannot be easily modeled as it is the case for errors presented in section 1.2.1. Hence, they remain as the main uncompensated error source. Their interference is usually presented as one of the main causes of precision degradation in the navigation solution. The purpose of this PhD is to study different approaches to try to minimize or mitigate the errors introduced by urban canyon phenomena.

Urban navigation phenomena

Multipath interference

Multipath is one of the largest GPS error sources, especially in weak signal environments [START_REF] Salytcheva | Medium accuracy INS/GPS integration in various GPS environments[END_REF]. In urban areas, the signal emitted from a satellite is very likely to get reflected and to follow different paths before arriving to the receiver. High buildings and tree foliage are typical reflecting surfaces.

Multipath signals are always delayed, since they travel longer distances than direct (LOS) satellite signals. They can cause the measured range to be too large or too small with respect to the true range, depending on the phase of the reflected signal.

The received signal will be thus composed by a sum of different and attenuated delayed replicas.

For satellite navigation purposes, only the direct signal is useful while the multipath components are considered as undesired signals. The expression for the nominal received composite signal with M i reflected components is obtained from (1.3) by setting i to the desired tracking satellite (i.e.

neglecting the presence of secondary satellites). It is very difficult to analyze the bias introduced by several simultaneous multipath signals on the estimation of the final TOA of the signal. The introduced error will depend on the power, phase and delay of each of the reflected signals [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF].

Considering a single multipath model and neglecting noise, the composite correlation function at the output of the receiver correlators can be expressed as

r i (τ ) = A 0 i 2 d i e -j[π∆f 0 i T b +(θ 0 i -θi )] sinc(∆f i T b )R x (τ -τ 0 i )+ A 1 i 2 d i e -j[π∆f 1 i T b +(θ 1 i -θi )] sinc(∆f 1 i T b )R x (τ -τ 1 i ), (1.15) where • (A 0 i , A 1 i
) represent the amplitudes of the corresponding LOS and multipath signal for the ith satellite,

• (∆f 0 i , ∆f 1 i ) represent the difference between the estimated and real Doppler frequencies for the LOS and multipath signal respectively,

• (θ 0 i , θ 1 i ) represent the phase of the corresponding LOS and multipath signal for the ith satellite,

• θi is the phase of the locally generated code replica for the ith satellite,

• sinc(•) denotes the cardinal sine function,

• R x (•) is the correlation function associated to the C/A spreading code,

• (τ 0 , τ 1 ) represent the time taken by the corresponding ith LOS and multipath signal to propagate from the satellite to the receiver.

The expression in (1.15) is formed by the resulting correlation between a local generated code replica with phase θ and both the direct signal and multipath delayed code. As a consequence the nominal correlation function to be used by the code discriminator (neglecting the phase and frequency error) will be deformed as depicted in Fig. 1.2. This deformation will entail a bias in the receiver computation of the direct signal parameters such as the time of arrival τ 0 , its amplitude A 0 and phase θ 0 . Hence, the posterior calculation of the position will be subjected to errors. In theory, pseudorange multipath error can reach magnitudes of about 0.5 of a code chip [START_REF] Parkinson | Global Positioning System: theory and applications[END_REF], i.e. 150m in the C/A case, depending on the receiver correlation technology.

For land applications, especially in kinematic mode, the magnitude of the multipath error quickly changes in a manner that is difficult to anticipate or to model mathematically [START_REF] Cannon | Analysis of a high performance C/A code GPS receiver in kinematic mode[END_REF]. In highly urbanized areas, LOS signals can be very weak with respect to the strength of the multipath signal or they can be completely blocked [START_REF] Mezentsev | Sensor aiding of HSGPS pedestrian navigation[END_REF]. In either situations, the receiver will track the reflected signal and a highly biased position estimation can be obtained. The magnitude of the induced error will no longer by limited to half a chip code (150m), because now the multipath is the main signal being tracked.

CHAPTER 1. GPS SATELLITE NAVIGATION

Cross correlation errors

Cross correlations effects are caused by the limited dynamic range of GPS Gold Codes and present a significant problem in the observation of weak GPS signals [START_REF] Glennon | A review of GPS cross correlation mitigation techniques[END_REF]. The use of 1,023 chip Gold-codes for the GPS C/A spreading-codes represents a compromise between the need for rapid acquisition and the cross correlation dynamic range of the spreading codes. Although Gold codes generally work well, there exist certain situations for which the tracking of weak GPS signals in the presence of other strong GPS signals can lead to cross-correlation errors. For these cases the use of these Gold codes represents a significant difficulty. A relatively short 1,023 chip Gold code with a chipping rate of 1.023 MHz and period of 1 ms was selected in order to minimize the code acquisition time, although a consequence of this choice is a cross correlation dynamic range of approximately 23 dB, which is not sufficiently large for some scenarios. Considering satellite i as the satellite being processed within the chosen channel and neglecting multipath and noise influence, the signal at the correlator output can be written as [START_REF] Van Dierendonck | Determination of C/A code self interference using cross-correlation simulations and receiver bench tests[END_REF],

r(τ ) = A i 2 d i e -j[π∆f i t+(θ i -θi )] R x (τ -τ i ) + ny k=1,k =i A k 2 d k T b e -j[π∆f k t+(θ k -θi )] c i (t)c k (t -τ k )dt, (1.16)
where

• the expression T b c i (t)c k (t -τ k )
dt represents the nominal cross correlation function between the desired satellite i and the secondary satellite k,

• ∆f k denotes the difference between the estimated ith Doppler frequency and the actual Doppler frequency of the kth satellite,

• θ k represents the kth signal phase.

Equation (1.16) shows that the cross correlation between satellite i and each other satellite k depends not only on their spreading-codes, but also on the Doppler frequency error. Indeed, if we considered that Doppler frequency f d,i for ith satellite is accurately estimated by the carrier tracking loops, the second term in (1.16) will also be determined by the relative Doppler-frequency difference between the ith and kth satellite [START_REF] Glennon | A review of GPS cross correlation mitigation techniques[END_REF]. In the case of GPS, the Doppler difference can result in several carrier cycles occurring within the integration time T b . As a result, prediction of any C/A code cross correlation must take into account the carrier phase difference between the codes. When the Doppler difference between satellite i and k is an integer multiple of 1 kHz, this mixing process will be relatively constant between successive C/A epochs and the cross correlated outputs will appear to be rotating at the same frequency as satellite i. Hence although the cross correlation samples rotate at the same frequency as the strong kth signal from which they are sourced, the frequencies are aliased into the bandwidth of satellite i. Data bit changes within the kth satellite signal can also cause a problem because since the codes are asynchronous, a data bit transition will typically occur within the code sequence for satellite k. This means that data-bits and relative Doppler effectively randomize the ideal cross correlations making them difficult to be predicted.

Low C/N 0
In urban scenarios, different types of obstacles, such as tree foliage or glass windows, can attenuate the received GPS signal causing a significant drop in the carrier-to-noise density ratio (C/N 0 ).The C/N 0 is the most fundamental parameter describing navigation signal quality. C/N 0 is an instantaneous measure of the ratio of the carrier power with respect to noise power density measured per Hertz of bandwidth (and similarly for dBW/ Hz). With a nominal noise floor spectral density of about -204 dB-Hz and minimum guaranteed signal power at -160 dB-W/Hz, the nominal C/N 0 magnitude is around 44 dBHz [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF]. When the C/N 0 goes beyond 28 dBHz approx., the signal is weak, likely erroneous, and highly susceptible to loss of lock by the receiver. This situation presents an important constraint for urban navigation where the quality of the received signal can be very low. In [START_REF] Macgougan | High sensitivity GPS performance analysis in degraded signal environments[END_REF] it was shown that pseudorange noise could increase by up to 1025m for signals degraded by 20 -25dB.

Carrier lock loops are specially susceptible to noise measurements. The amount of tolerable frequency error during the total dwell time depends on the length of coherent integration and the type of carrier tracking performed [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF]. A common approach to deal with these types of situations is to use high sensitivity techniques that enable one to wider the C/N 0 range for which signals can be tracked. This solution especially adapted to degraded scenarios is presented in the next section.

Reduced visibility over the GPS satellite constellation

The lack of redundant measurements in urban scenarios presents a main navigation problem. High buildings block the visibility of medium/low elevation satellites and tunnels completely disrupt the GPS based navigation. In this way, partial (less than 4 visible satellites) or total (no visible satellites) outages prevent the continuity over the position estimation to be guaranteed. This is a critical issue for applications such as safety of life.

High Sensitivity GPS

The accuracy degradation caused by signal masking (low C/N 0 ) and poor geometry can be mitigated to some extent by the use of innovative High Sensitivity (HS) GPS technology [START_REF] Peterson | Measuring GPS signals indoors[END_REF]. Unlike conventional GPS receivers, High Sensitivity GPS (HSGPS) receivers make measurements in signal conditions where conventional sensitivity receivers fail. HS receivers are capable of tracking and acquiring signals in challenging environments, which often include the interior of buildings, under heavy foliage and in urban canyons [START_REF] Chansarkar | Acquisition of GPS signals at very low signal to noise ratio[END_REF]. HS receivers are specially designed to track low power GPS signals. The L1 C/A-code repeats every millisecond. This can be used by the GPS receiver, so that the signal can be integrated for extended periods of time in order to obtain a higher Signal-to-Noise Ratio. The longer integration time, the better the reduction of the noise power at the output of the correlator. The increase in the coherent/non-coherent integration time allows signals with power as low as -188dBW to be processed. However, this technique is still affected by multipath, jamming or cross-correlation errors. Longer coherent integration is possible, if the navigation bits are known a priori, but this process is still restricted by the residual errors due to receiver and satellite movements during the integration interval [START_REF] Salytcheva | Medium accuracy INS/GPS integration in various GPS environments[END_REF]. In fact, here relies one of the main HSGPS limitation: under very low signal-to-noise ratios, the bit error rate associated to the demodulation of the navigation message is extremely high. hence, without indispensable information such as ephemerids or the satellite clock correction, the navigation solution cannot be computed.

The Assisted GPS technology [START_REF] Kubrak | Hybridisation of a GPS receiver with low-cost sensors for personal positioning in urban environments[END_REF] appeared in the last years as a solution to enhance HSGPS performances and its tracking of weak signals. Navigation data is transmitted through the mobile network to the receiver. In this way, a navigation solution can be provided even in very degraded scenarios where the HSGPS technique fails. However, this technique is out of the scope of this thesis.

Conclusion

The general structure of the Global Positioning system (GPS) was introduced in this chapter. The transmitted signal properties were presented. Possible measurements that can be obtained at the receiver output were described together with their error sources. Special attention was given to urban canyon phenomena were the received signals are susceptible to multipath interference, crosscorrelation effects and/or signal attenuation. All these concepts will be of a major importance in the next chapters. Indeed, the scope of this thesis is to propose different approaches to deal with the lack of integrity presented by the erroneous signals (hereafter referred as outliers) in urban environments.

Chapter 2 Gyroscopes, accelerometers and wheel speed sensors (WSS) are among the most common DR sensors used for land vehicle navigation. Other sensors such as magnetometers or barometers might also be included in the navigation system to enhance performances. The Inertial Navigation System (INS) [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF] is an especially well known example of dead reckoning navigation combining accelerometers and gyroscopes. If WSS information is available for two opposite wheels, a WSS-based navigation system is also possible [START_REF] Zhao | Vehicle location and navigation systems[END_REF].

Dead reckoning navigation

In many aspects, dead reckoning systems are a perfect complement to GPS. They are self-contained and not affected by atmospheric propagation errors. However they need an absolute navigation system such as GPS to provide an initial position estimation and to control their time diverging errors. DR systems are particularly interesting in urban scenarios where GPS availability and reliability might be frequently jeopardized. DR sensors provide an additional source of measurements in case of partial or total GPS outage, and can help to mitigate or filter outlier presence.

This chapter begins with a presentation of the dead reckoning sensors to be used within this thesis.

As one of the main challenges of DR navigation involves coordinate frame transformations between local and absolute reference frames, section 2.5 gives a brief insight into the different navigation frames to be studied herein. Equations of motion are given in section 2.6 and 2.7 for the INS and WSS-based navigation system respectively. Based on the error sources associated to each DR navigation system, their propagation error equations are also specified.

Inertial Sensors

Gyroscope

A gyroscope (also referred as gyro) is an instrument used to measure the rate of rotation or integrated heading change of a platform. Quite a few different types of gyros are available, ranging greatly in price and stability. Details for the different existing gyros can be found in for example in [START_REF] Grewal | Global Positioning Systems, Inertial Navigation and Integration[END_REF] and [START_REF] Kayton | Avionics Navigation systems[END_REF]. Gyros are classified into gimbaled or strap-down varieties, with gimbaled gyros maintaining a fixed orientation in an inertial frame (relative to distant galaxies). Strap-down gyros measure rotation on a fixed plane with respect to the vehicle, which is not generally on the plane orthogonal to the gravity vector. Low cost sensors, as the one being used within this thesis, generally belong to this last category.

Sensors are affected by many errors that vary with time, temperature, as well as with the motion experienced by the sensor. The magnitude of these errors will depend on the quality of the used sensors. The theoretical output model for the gyros can be expressed as [START_REF] Jekeli | Inertial navigation systems with geodetic applications[END_REF]:

ω out = (1 + SF g )ω real + δSF g ω real + b ω + δT p g + s mag + m g alig + n g , (2.1)
where

• ω out is the rotation rate delivered by the gyro,

• SF g is the scale factor affecting the true rotation rate,

• ω real is the true rotation rate,

• b ω is the bias affecting the gyro measurement,

• δSF g is the gyro scale factor error,

• δT p g is the gyro temperature sensitivity effect,

• s mag is the magnetic sensitivity term,

• m g alig is the misalignment error

• n g is the gyro noise with n g ∼ N (0, σ 2 g ).

Accelerometers

An accelerometer measures platform acceleration, which can be integrated to give velocity, and double integrated to give distance traveled. Accelerometers are generally based on observing the displacement of a suspended mass caused by inertia [START_REF] Stephen | Development of a multi-sensor GNSS based vehicle navigation system[END_REF]. This mass can be classically suspended by either a horizontal spring or vertical pendulum. Once a force is applied to this mass, the movement of the mass would be sensed through capacitance, an optical method, or usually, by the use of a viscous fluid medium. If the measured force is divided by the mass of the suspended mass, the final specific force is obtained. This specific force f is the true acceleration measure by the accelerometer. Similarly to the model used for the gyroscope in (2.1), the accelerometer output model is typically

given by [Jek00]:

f out = (1 + SF a )f real + δSF a f real + b f + δT p a + e + m a alig + n a , (2.2)
where

• f out is the specific force delivered by the gyro,

• SF a is the scale factor affecting the force,

• f real is the true specific force,

• δSF f is the accelerometer scale factor error,

• b f is the bias affecting the accelerometer measurement,

• δT p a is the accelerometer temperature sensitivity effect,

• e is the error due to loss of elasticity,

• m a alig is the misalignment error,

• n a is the accelerometer noise with n a ∼ N (0, σ 2 a ).

Adopted noise models

Modeling and estimating all the errors associated to both the accelerometer and gyro presents a very high cost. First of all, a large number of measurements coming from a complementary system such as GPS would be necessary in order to guarantee the observability of all the parameters [START_REF] Petovello | Real-Time Integration of a Tactical-Grade IMU and GPS for High-Accuracy Positioning and Navigation[END_REF]. This condition could be hardly fulfilled in urban scenarios. The estimation of such a large number of errors would jeopardize the stability and real-time response of the system. Moreover, while working with low cost sensors, most of the errors would be indistinguishable from the noise component. Therefore, the sensor model to be used hereafter will only be affected by a bias and noise term. The simplified measurement model for the gyroscope and accelerometer is thus expressed as [START_REF] Giremus | Apports des techniques de filtrage non linéaire pour la navigation avec les systèmes de navigation inertiels et le GPS[END_REF]:

ω out ≈ (1 + SF g )ω real + b g + n g ← gyro model (2.
3)

f out ≈ (1 + SF a )f real + b a + n a ← accelerometer model (2.4)
where any initial run-to-run bias for the gyro is considered to be compensated during the sensors calibration stage. A first order Markov model is used for the bias characterization as follows:

ḃg = - 1 T g b g + w bg , (2.5
)

ḃa = - 1 T a b a + w ba , (2.6)
where

• b a and b g are not necessarily equal to b f and b ω respectively, but they represent all the possible bias sources affecting the accelerometer and gyro measurements,

• T a and T g correspond to the process correlation time,

• (w ba , w bg ) are gaussian noises described by w ba ∼ N (0, σ 2 ba ) and w bg ∼ N (0, σ 2 bg ) respectively. These parameters account for the temporal variability of the time-varying parameters used in the theoretical model [(2.1),(2.2)] such as temperature sensitivity effect, scale factors, etc. Furthermore, considering that the time constants T a and T g are usually rather large, the biases can be approximated by ḃa ≈ w ba and ḃg ≈ w bg .

A small parenthesis is done on the description of DR sensor models, to introduce a tool that will be employed in the following chapters to obtain an accurate characterization of the accelerometer and gyroscope error terms. According to the sensors' output models given in (2.3)-(2.4) and (2.5)-(2.6), four noise terms need to be evaluated: (n g , w bg ) for the gyroscope and (n a , w ba ) for the accelerometer.

With this purpose, the Allan Variance technique will be briefly presented in the following section.

Allan Variance

The Allan variance [All66] is a method of representing the root mean square (RMS) random-drift error as a function of averaging time. It is easy to compute and relatively simple to interpret and understand. The Allan variance method can be used to determine the characteristics of the underlying random processes that give rise to data noise. In particular it can be used to characterize the different types of noises affecting the inertial-sensor data as shown in [START_REF] El-Sheimy | Analysis and modeling of inertial sensors using Allan variance[END_REF] and [START_REF] Zhang | Allan variance analysis on error characters of MEMS inertial sensors for an FPGA-based GPS/INS system[END_REF]. To compute the Allan variance, the output of the idle sensor is measured. The idea is that the data set is divided in smaller data sub-sets (or sample groups) and the mean of each of these groups is computed. By comparing the means of consecutive sample groups, the variance associated to the whole data set can be calculated. The obtained variance will necessarily depend on the length, or averaging time, of the sub-sets. Consider for instance that fig. 2.2 represents the output of an idle sensor.

σ 2 AL (α) = 1 2 bn+1 -bn 2 , (2.7)
where denotes an infinite time average which, in practice, is necessarily approximated by a finite sum including the N data samples. It is shown in [START_REF] Tehrani | Ring laser gyro data analysis with cluster sampling technique[END_REF] that by varying the integration time α of the sample groups, different random processes can be detected. A direct relation can be established between the Allan variance and the power spectrum density of the desired random process. In particular, the Allan deviation characteristic curve for the sensor white noise n (where n stands for n g or n a ) can be derived by a log-log calculation:

log(σ n AL (α)) = - 1 2 log(α) + log (σ AL (1)) , (2.8)
while for the noise w b driving the first order Markov model associated to the sensor bias (where w b stands for w bg or w ba ) the expression is the following:

log(σ b AL (α)) = 1 2 log(α) + log σ AL (3) (3)
.

(2.9)

A typical curve for the Allan variance as a function of the sample group length α is given, in logarithmic scale, in fig. 2.3. It can be easily deduced from (2.8) that if a tangent is drawn for the logarithmic Allan variance curve with slope -1/2 and its value in α = 1 is evaluated, the white noise n spectral density is obtained (and straightforwardly the noise variance) [START_REF] Zhang | Allan variance analysis on error characters of MEMS inertial sensors for an FPGA-based GPS/INS system[END_REF]. Similarly, from (2.9), if a tangent is drawn for the Allan variance curve with slope 1/2 and its value in α = 3 is evaluated, the white noise spectral density for the w b noise driving the Markov model random walk is obtained. This method will be used in the next chapter to obtain an accurate characterization of the sensor error terms.

Odometers/Wheel Speed Sensor

An odomoter is a sensor that measures the distance traveled by a vehicle, or possibly by the individual tires, with respect to an initial position [START_REF] Zhao | Vehicle location and navigation systems[END_REF]. A wheel odometer or wheel speed sensor (WSS)

typically measures the number of rotation counts (pulses) generated by a rotating wheel. This mechanical motion (i.e. the rotation counts) is then converted into an electrical signal whose frequency is proportional to the rotational velocity of the wheel.

The Hall-effect odometer is one of the most common WSS installed on vehicles. It allows the sensing of low wheel speeds better than many other simple mechanical or pulse detection sensors.

The Hall effect occurs when a charge carrier moving (i.e a current) through a material experiences a deflection because of an applied magnetic field. This deflection results in a measurable potential difference across the side of the material that is transverse to the magnetic field and the current direction. wheel radius, the velocity of the vehicle (in m/s) and therefore the distance traveled can be computed.

ODOMETERS/WHEEL SPEED SENSOR

Errors associated to the speed estimation obtained from the WSS will be mainly related to slippage effects and to errors in the nominal tire circumference as explained in [START_REF] Stephen | Development of a multi-sensor GNSS based vehicle navigation system[END_REF] and [START_REF] Svenzén | Real time implementation of MAP aided positioning using a Bayesian approach[END_REF]. Slippage occurs if the vehicle accelerates or decelerates too rapidly or travels on a snowy/icy/wet surface. The wheel diameter may vary due to changes in the tire pressure over time, wear of the tire or changes in temperature. At the same time, wheel rotation causes a centrifugal force on the tire material that tends to stretch it out as shown in Fig. 2.5.

Considering the different above-mentioned types of errors affecting the wheel velocity determination as well as the noise associated to the specific sensing of the rotational velocity, the output model to be used for the WSS is expressed as:

v wss = ω w (R w + δR w ) + n w (2.10)
where

• v wss is the velocity obtained from the WSS in m/s,

• ω w is the wheel rotational rate in rad/s,

• R w is the wheel nominal radius,

• δR w represents a variation on the nominal radius due to temperature, pressure, speed...etc,

• n w describes both the noise associated to the senor and possible errors introduced by the wheel slippage.

The radius error is described as white noise given by:

δR w = n r , (2.11)
where n r ∼ N (0, σ 2 r ). It has to be noted that except from extraordinary situations such as flat tire, the variation in the wheel radius is really small. For instance, modern tires have fairly constant size and shape and they are designed to dissipate large amounts of heat caused by road friction. At the same time, specially adapted tires for snowy or iced conditions and the ESP (Electronic Stabilization Program) vehicle option help to reduce the slippage effect. In practice, the slippage effect can be almost neglected for civilian land vehicle navigation (see [START_REF] Enquist | Aspects of high precision estimation of vehicle dynamics[END_REF] for more details). Most vehicles are nowadays equipped with WSS to track the total distance traveled by the car.

WSS are used by the anti-lock braking system (ABS) and traction control system. ABS WSS data will be exploited in this thesis. The ABS is becoming a standard feature in most vehicles, providing easily accessible odometric data at no additional cost. Indeed, messages containing the wheels' angular rate velocity will be read from the ABS control area network (CAN) bus. Every car equipped with an ABS contains an "on board diagnostic" (OBD) port which provides the necessary interface between the CAN bus and the user. However, standards followed by the data transfer through this interface are not directly compatible with personal computers (PCs). An "interpreter" device (Fig. 2.6) is needed to act as a bridge between the OBD port and the nominal serial RS232 PC port. Most ABSs provide two to four-wheel control. These different measurements can be exploited to compute speed, distance traveled and heading by a technique called differential odometry. This navigation strategy is presented in section 2.7.

Coordinate frames

All the above mentioned sensors are usually rigidly mounted on the vehicle so their measurements are given in the vehicle-dependent or mobile frame. A frame transformation must be done to the navigation or ECEF frame to enable data to be exploited by the user as positioning information. In fact, this aspect represents one of the main difficulties of the dead reckoning navigation system.

Before discussing the position calculation, the different coordinate frames are going to be presented.

Inertial Frame (i-frame)

An inertial frame is considered as non-rotating and non-accelerating frame relative to distant stars.

The Earth's centre of mass is its origin, with X i pointing toward the mean vernal equinox and Z i parallel to the earth's instantaneous spinning axes (pointing toward the north pole). Y i is aligned to complete a right-handed frame. The frame coordinates are depicted in Fig. 2.7(a) as i = (X i , Y i , Z i ).

Earth-Centered Earth-fixed frame (ECEF or e-frame)

The ECEF, described by its (X e , Y e , Z e ) axes, is linked to the Earth. It is deduced from the inertial reference frame by applying a rotation about the Z i = Z e axis at the earth's rotation rate ω ie . The Z i = Z e assumption is done by neglecting effects mainly due to polar motion [START_REF]Navstar GPS space segment/navigation user interfaces[END_REF]. The X e axis is given by the intersection between the equatorial plane and the plane containing the Greenwich meridian.

Navigation frame (n-frame)

The Navigation frame n = (n, e, d) has its origin in the vehicle's centre of mass. The north axis n points toward the true north, the east axis e points toward the east and the down axis completes the right handed coordinate system pointing toward the interior of the Earth, but not necessary to the earth centre. It can be obtained by two consecutive rotations around the X e and Z e axes of the e-frame by magnitudes of geodetic latitude, λ, and longitude, ϕ. The navigation frame is presented in Fig. 2.7(b).

Mobile frame(m-frame)

The mobile reference frame is attached to the mobile. In strap-down DR systems, as the ones used

Inertial Navigation system

Inertial navigation is an autonomous process of computing position by doubly integrating the acceleration of a point, whose position is to be determined [START_REF] Jekeli | Inertial navigation systems with geodetic applications[END_REF]. It was developed for the second world war and mainly conceived for aerial and spatial applications. This system relies on the measurements provided by sensors contained in an Inertial Measurement Unit (IMU) [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF]. This IMU is generally composed of three accelerometers and gyroscopes.

There are two classifications of inertial systems: gimbaled and strapdown. A gimbaled INS deals with the physical realization of the navigation frame using a free axis gyrostabilizer platform with orthogonally placed accelerometers. Gimbaled gyroscopes are used to aligned the platform frame with a specific navigation coordinate system. The isolation of the inertial platform from rotations of the host vehicle can be exploited to eliminate many sensor error sources and to achieve very high system accuracy [START_REF] Titterton | Strapdown Inertial Navigation Technology. Radar, sonar and navigation[END_REF]. However, gimbals are very sophisticated electromechanical assemblies that are delicate and expensive to manufacture. For low cost land applications, as the one studied within this thesis, strap-down systems are normally used. In strap-down systems, the sensors are fixed to the vehicle itself and not to a stabilised platform. Strapdown gyroscopes are not used to keep the accelerometer input axes stabilized, but they are used to maintain a coordinate transformation between the vehicle and navigation frame. Inertial sensors for strap-down systems experience much higher rotation rates than gimbaled systems and, consequently, they deliver poorer accuracy than their gimbaled counterparts.

An INS is composed of an IMU and a computer [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF]. The IMU provides the raw measurements at the output of the inertial sensors. These measurements, including rotation rates and specific force measurements, are converted by the computer into position, velocity and attitude information.

A 2D INS is particularly studied in this thesis. The IMU is composed by a 1 axis gyroscope measuring the yaw rate and by 2 axes accelerometers in the along track and cross track direction.

Altitude h is considered almost constant and known. All expressions given hereafter will concern a 2D

INS. In fact, the choice for a reduce inertial navigation system was motivated by the low accuracy of the available sensors. For sensors aligned with low dynamic axes, such as the down accelerometer or the gyros measuring rotational rates around the (x m , y m ) mobile coordinates, no useful information could be generally extracted from the noisy output measurements.

Notation

In order to fully understand the INS equations presented in the next section, the following notations relating parameters and their associated coordinate frames need to be introduced:

• x a is the projection of vector x into frame a.

• ω a bc is the rotation rate vector of frame c, relative to frame b, expressed in the a frame.

• R a2b is the rotation matrix from the a to b frame.

Equations of motion

Once the angular rates and accelerations are obtained from the IMU, the calculation algorithm can be executed. Firstly, having the initial navigation information with respect to the navigation frame, the angular rates are integrated to obtain the new orientation of the IMU. Secondly, using this information, the accelerations are rotated into the n-frame, where they are twice integrated to obtain velocity and position increments. The frame transformation relating velocity and position in the navigation frame is given in (2.12). The geodetic position [λ, ϕ] considering the velocities in the navigation frame v n = [v n , v e ] are known is obtained according to:

ṗn =   λ φ   =   1 R λ 0 0 1 Rϕ cos(λ)     v n v e   , (2.12)
where

• ȧ denotes the rate of a in a continuous time domain,

• λ is the latitude rate,

• φ is the longitude rate,

• R λ is the earth radius of curvature in a meridian at a given latitude,

• R ϕ is the transverse radius (considering the WGS84 model for which the earth is an ellipsoid).

The Earth's radii are obtained as follows

R λ = a(1 -e 2 ) 3 2 1 -e 2 sin 2 (λ) , R ϕ = a 1 -e 2 sin 2 (λ) , (2.13)
where a is the semi-major axis and e the eccentricity of the WGS-84 ellipsoid. Considering a two dimension navigation solution, the velocity in the navigation frame is obtained as

vn = R m2n f a + a e + a c , (2.14)
where

• vn is the rate of the velocity vector in the navigation frame,

• R m2n is the frame transformation matrix from the mobile to the navigation frame,

• f a = [f u , f v ]
T is the 2D specific force vector in the mobile frame (i.e. the accelerometer output),

• R m2n f a = [f n , f e ] T represents the specific forces projected in the 2D navigation north-east frame,

• a e = -ω n ei ∧ ω n ei ∧ v n are the centrifugal accelerations induced by the Earth rotation, also called gravity anomaly parameters. The operator ∧ denotes the cross-product function. This term will be neglected in the following because it is generally masked by the uncertainty of the measurement itself.

• a c = (ω n en + 2ω n ie ) ∧ v n represents the Coriolis acceleration (to be considered in 2D)

Neglecting the gravity anomaly parameters introduced by the centrifugal forces, the developed expression of eq. (2.14) can be written in a vector form as:

vn =   vn ve   = R m2n f a +   -v e ( φ + 2ω ie ) sin(λ) v n ( φ + 2ω ie ) sin(λ)   , (2.15)
where the last left-size term represents the Coriolis acceleration and ω ie is the Earth rotation rate (≈ 7.2924 • 10 -5 rad/s = 4.17 • 10 -3 deg/s).

The transformation frame matrix R m2n is computed from the yaw rate provided by the gyro such as:

R m2n =   cos(ψ) -sin(ψ) sin(ψ) cos(ψ)   (2.16)
where ψ represents the yaw rate measured by the along track gyro (i.e the gyro aligned with the X m coordinate of the mobile frame). In theory, the yaw rate angle ψ is not directly obtained from the gyro output but according to the following expression:

ψ = ω m nm = ω m im -R n2m (ω n ie + ω n en ), (2.17)
where ω m im is gyro output, and ω n ie is the inertial rotation rate of the Earth expressed in the navigation frame so

ω n ie = w ie [cos(λ), 0, -sin(λ)] T . (2.18)
The maximum rate for eq. ( 2.18) is obtained either when the vehicle is at the equator or at the Earth's pole. The transport rate ω n en of the navigation frame relative to the Earth is application dependent and it is given by

ω n en = [ φ cos(λ), -λ, -φ sin(λ)] T . (2.19)
However, in practice, the second term in the right hand side of equation (2.17) can be considered negligible when compared with typical noise values affecting low cost gyros used for land vehicle navigation. It was shown in [START_REF] Kubrak | Hybridisation of a GPS receiver with low-cost sensors for personal positioning in urban environments[END_REF], that this term can be considered as acting as a constant bias of magnitude 5.5×10 -3 deg/s, which is generally far below the noise affecting the gyro output. Therefore, the rotation rate of the mobile frame with respect to the navigation frame expressed in the mobile frame ω m nm is approximated to the gyro output as:

ψ = ω m nm ≈ ω m im . (2.20)
The yaw angle is thus updated at time t as

ψ(t) = ψ(t -1) + ω m im ∆t (2.21)
where ∆t is the sampling period.

INS mechanization

The mechanization of the inertial navigation system consists of the different steps followed by the inertial measurement unit in order to provide a navigation solution. These steps, for a 2D solution, are mainly summarized by the differential equations given in (2.12) and (2.15), and can be separated in four main categories:

1. Attitude update

Error equations

No error considerations were included in the theoretical motion models presented in the previous section. However, it is generally critical to have an estimate of these errors for the otherwise time diverging INS position solution. Assuming that the vehicle position is accurately known during the initialization process, the precision of the computed solution will mainly depend on the quality of the inertial sensors. As already presented in eqs. (2.3) and (2.4) the gyro and accelerometer outputs are considered to be affected by a white noise and a time varying bias that summarizes the contribution of different error sources to the sensor measurements. Gyro errors will impact the Rm2n matrix computation by introducing frame transformation errors. Accelerometer errors will be consequently integrated to compute the final position. Errors introduced by the sensor biases on the velocity or position computation can be derived from eq. ( 2.12) and (2.15). It is important to notice that the successive integrations needed to compute velocity and position result in an approximative linear and parabolic growth of the errors associated to the sensor biases. Just to give a clarifying example of the impact the accelerometer and gyro biases have on the velocity and position calculation lets consider the following ultra-simplified error equations (neglecting the Earth rotation effect and the dependence between the north and east axes):

     δ ṗne δ vn δ ψ      Xs =      0 I 2 0 0 0 -f n 0 1 R 0      F      δp ne δv n δψ      +      0 b n a b g      B (2.22)
where

• δ• denotes the difference between the actual value and the calculated (or measured) quantity,

• δp ne denotes the position errors in the north-east frame so δp n = Rδλ and δp e = R cos(λ)δϕ,

• δv n and δψ denote the corresponding velocity and angle errors,

• I 2 is a 2 × 2 identity matrix,
• R accounts for an average Earth radius,

• -f n =   -f e f n 
 is the vector including the accelerometer forces in the navigation frame,

• b n a represents the accelerometer bias vector in the navigation frame,

• b g accounts for gyroscope bias.

Considering null initial conditions the analytical solution is given by [START_REF] Giremus | Apports des techniques de filtrage non linéaire pour la navigation avec les systèmes de navigation inertiels et le GPS[END_REF]:

X s (t) = t t 0 exp( F (t -u)) Bdu (2.23) so, δp ne (t) ≈ b n a ∆T 2 2 + f n b g ∆T 3 3 , (2.24) δv(t) ≈ b n a ∆T + f n b g ∆T 2 2 , (2.25)
where ∆T = t -t 0 . Generally speaking, for short periods of time, the inertial navigation errors (2.24)-(2.25) appear as exponentially divergent. It is thus crucial to estimate the biases affecting the INS.

Linear equations are usually employed to describe the INS error dynamics. Therefore, a first order approximation of the INS differential navigation equations is done. Given the real complexity of the system, no full derivation is given here. The adopted 2D model for INS error state equations used in this thesis is directly derived from the 3D model given in [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF]. The linearized state transition model for position, velocity, heading angle and sensor biases is summarized as:

           δ ṗn δ vn δ ψ ḃa ḃg            IN S =            F pp F pv F pϕ 0 2×2 0 2×1 F vp F vv F vϕ F vba 0 2×1 F ϕp F ϕv F ϕϕ 0 2×2 F vbg 0 2×2 0 2×2 0 2×1 -diag(1/T a ) 0 2×1 0 1×2 0 1×2 0 0 1×2 -diag(1/T b )            F IN S t            δp n δv n δψ b a b g            IN S +            0 2×2 0 2×1 0 2×2 0 2×1 D a R m2n 0 2×1 0 2×2 0 2×1 0 1×2 σ g 0 1×2 0 0 2×2 0 2×1 D ba I 2×2 0 2×1 0 1×2 0 0 1×2 σ bg            BINS t w, (2.26)
where

• δp n = (δλ, δϕ) denotes the INS position errors in the navigation frame,

• δv n = (δv n , δv e ) denotes the INS velocity errors in the navigation frame,

• δψ is the heading (or yaw) angle error,

• 0 i×j denotes a i × j zero matrix,

• D a =   σ x a 0 0 σ y a 
 is a diagonal matrix containing the noise standard deviations (stds) for the accelerometers aligned with the x and y axis of the mobile frame,

• D ba =   σ x ba 0 0 σ y ba 
 is a diagonal matrix containing the Markov model noise stds for the accelerometers aligned with the x and y axis of the mobile frame,

• w is a gaussian noise vector so w ∼ N (0, I 6×1 ).

A detailed expression for all sub-matrices is found in [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF]. Please note that the above presented sub-matrices are a reduced version of the 3D sub-matrices given in [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF], where just the 2D pa-

Navigation based on Wheel Speed Sensors

This section describes the main elements of differential odometry [START_REF] Stephen | Development of a multi-sensor GNSS based vehicle navigation system[END_REF]. Differential odometry is a technique to provide both distance traveled and heading rate information by integrating the outputs from two odometers, corresponding to a pair of front or rear wheels. When a vehicle turns, the differential allows the inside tire to travel a shorter distance than the outside tire without greatly increased slippage. The individual wheel speeds would therefore vary from the actual along track speed while the vehicle is turning. The vehicle velocity is computed by averaging the two wheel velocities. The heading change is computed as the difference between the wheel' velocities divided by the axle length (i.e. length between wheels). A 2D position solution can be then calculated from the average speed and heading measurements. Figure 2.10 shows WSSs located on the front and rear wheels. The first index of the different variables refers to the front f or rear r axes whereas the second index corresponds to the left l and right r sides of the car. Consequently, the wheel radii (resp. angular velocities) are denoted as R rl , R rr , R f l and R f r (resp. ω rl , ω rr , ω f l and ω f r ). The other notations used in Fig. 2.10 are L for the length between wheels and ψ for the vehicle yaw rate.

This work focuses on velocity and yaw rate calculations using rear wheels. This choice is motivated by the fact that acceleration and deceleration have less effect on the output of the non driven wheels [START_REF] Zhao | Vehicle location and navigation systems[END_REF].

Assume first that the wheel radii are constant and known. The mean speed of the vehicle can be computed as [START_REF] Svenzén | Real time implementation of MAP aided positioning using a Bayesian approach[END_REF]:

V = ω rl R rl + ω rr R rr 2 + n V . (2.27)
where n V ∼ N (0, σ 2 V ) represents the noise due to errors in the sensing of the angular velocities,

ψ = ω rl R rl -ω rr R rr L + n ψ . (2.28)
where n ψ ∼ N (0, σ 2 ψ ) considers the impact in the yaw angle calculation of the same error sources described for n V . Assuming a constant vehicle altitude and a correct initialization of the yaw angle, changes in the position can be computed as:

  v n v e   =   V cos(ψ) V sin(ψ)   ,
(2.29)

Error equations

Errors in the wheel radii will have a strong impact on localization accuracy when getting propagated through the angle and velocity expression. Indeed, any non corrected error will result in an accumulative increasing error. In particular, as the yaw rate is computed from the difference between the two wheel velocity, a non corrected wheel radius error, will entail a very high proportional error, so that when integrated over time, the final yaw angle will substantially differ from the really yaw angle.

By using (2.28), the resultant bias produced at time t by a radius error can be written as [START_REF] Carlson | Practical position and yaw rate estimation with GPS and differential wheelspeeds[END_REF]:

bias ψ ≈ t 0 ω rl t δR rl t -ω rr t δR rr t L dt, (2.30)
where δR rr t and δR rl t are modeling the radius variations in the time interval [0, t]. Therefore the interest of modeling and estimating the radii errors. The linearized error equations proposed for the differential odometry navigation, derived from eqs. ((2.27),(2.28),(2.12),) , are expressed as follows

              δ λ δ φ δ V δ ψ δ Ṙrl δ Ṙrr               W SS =               0 0 cos(ψ) R λ +h V -sin(ψ) R λ +h 0 0 V tan(λ) 2 Rϕ+h 0 sin(ψ) (Rϕ+h) cos(λ) V cos(ψ) Rϕ+h 0 0 0 0 0 0 ω rl 2 ω rr 2 0 0 0 0 ω rl L -ω rr L 0 0 0 0 0 0 0 0 0 0 0 0               F W SS t               δλ δϕ δV δψ δR rl δR rr               W SS +               0 0 0 0 0 0 0 0 σ V 0 0 0 0 σ ψ 0 0 0 0 σ rl 0 0 0 0 σ rr               BW SS t w, (2.31)
where

• δV is the WSS along track velocity error,

• δψ is the WSS yaw rate error,

• δR rl represents the error in the nominal radius of the rear-left wheel,

• δR rr represents the error in the nominal radius of the rear-right wheel,

• w is a gaussian noise vector so w ∼ N (0, I 4×1 ).

To overcome WSS accuracy limitations to the yaw angle computation, chapter 4 proposes a WSSbased navigation system augmented by a 1 axis gyro. In this way, the WSSs only provide speed measurements and the gyro provides the necessary yaw angle data. This is one of the most commonly used dead reckoning navigation approaches for land vehicle navigation. The only difference with the above presented model is that ψ is directly measured from the gyro (as in the INS) and not obtained as a function of the wheels' angular velocities. The gyro error model given in (2.5) is thus used to define the error equations.

Conclusions

This chapter addressed different DR-based navigation systems for 2D land vehicle navigation. The inertial navigation system (INS) was presented along with a Wheel Speed Sensor(WSS)-based navigation system. Motion and error equations were given for both systems. All the explained DR principles given herein will be used in the next chapter to define a filtering model to be applied to hybrid navigation solutions.

Chapter 3

Hybrid navigation systems: integrating

Reliability and continuity over the navigation solution cannot be generally ensured for the GPS.

Several phenomena, such as lost of visibility over the satellite constellation, referred as outages, or highly corrupted signals, referred as outliers, can seriously degrade the navigation performance.

Additional measurement sources are thus necessary to compensate for these effects. This chapter exploits information given by inertial sensors and on-board Wheel Speed Sensors (WSS) to propose hybrid navigation techniques.

The Extended Kalman Filter (EKF) is traditionally used within the navigation community. It has a low computational cost and it provides general acceptable performances. However, it can lead to suboptimal results if the first order linearization that must be applied to the system model entails large linearization errors. As a response to this limitation, an enhanced solution started gaining a wide acceptance during the last decade: the Unscented Kalman Filter (UKF). The properties of both Kalman-based filtering techniques will be explored in this chapter and its performances contrasted.

Integration strategies

There exists several strategies to integrate GPS and dead-reckoning data into a common system. They all depend on the measurements available at the GPS level and the way they are combined with dead-reckoning data. The best choice for an integration strategy is generally application-dependent and it is commonly determined by the quality of the sensor data, the navigation scenario and the affordable system complexity. The main principle shared by all the integration strategies is that GPS updates are utilized to calibrate the dead-reckoning sensors.

The most intuitive approach to a hybrid GPS/DR system would be to estimate the vehicle cinematic parameters. However, high vehicle dynamics would demand the estimates to be updated at a very high frequency and therefore a very short filter time response would be needed. To overcome this expensive filtering approach, DR navigation errors are usually estimated. As errors experienced lower dynamics than navigation parameters, the filtering conditions are relaxed. There are four main types of integration strategies: the loose, tight, ultra-tight and the deep integration [START_REF] Jekeli | Inertial navigation systems with geodetic applications[END_REF]. In the loose strategy (also referred as "loosely coupled") GPS position and velocity information are exploited, while the tight strategy (also referred as "tightly coupled") uses the raw GPS measurements, such as the pseudorange or Doppler measurements, for the integration approach. In the ultra-tight strategy, integration is done at a very high level. GPS updates are obtained at the correlator outputs and they are utilized to calibrate the DR system. At the same time, DR is used to aid the GPS receiver tracking loops during interference or other degraded signal conditions [START_REF] Sennott | Robustness of tightly coupled integrations for real-time centimeter GPS positioning[END_REF]. An extension of this concept is the deeply coupling method. In this case the integrating navigation filter is implemented as one element of the receiver tracking loop. Using inertial information in combination with the in-phase and quadra-phase signals of the receiver signal processing an optimal controlling value for the numerical oscillator can be computed. Optimisations of signal tracking performance are the goal of this method.

However, these last two types of integration methods require access to the receiver firmware, or at least to the tracking loop information. They are usually implemented by equipment manufacturers and thus they are not within the scope of this work. Therefore, further discussions will only concern loose and tight integration. Both strategies will be presented for an open-loop configuration.

Loose integration

The loose integration scheme fuses GPS and DR at the lowest level. It works with the position and velocity estimations at the output of the GPS navigation filter, and compares them to its analogous DR estimations. The GPS filter independently processes the tracked satellites and provides a position For this integration method GPS and DR position solution are computed in separated filters. This division enables a low complexity system implementation while precision is sacrificed. A minimum of four visible satellites is needed for the GPS navigation solution to be determined and the DR growing errors to be compensated. This condition cannot be guaranteed in urban scenarios or under situation where the received signal is degraded. Moreover, if an outlier is present, no fault detection techniques can be applied at the measurement level. Therefore an undetected error would get easily propagated to the position solution.

Tight integration

Tight integration works with raw GPS measurements (such as pseudorange or Doppler measurements) which are used as an input to the "coupling filter". In this strategy GPS measurement, and not navigation parameters, are independently combined with the DR data. Indeed, DR position cinematic parameters are used to estimate the received GPS raw measurements. The "coupling filter" estimates DR errors, and corrections are applied directly on the DR navigation solution. A tight integration scheme is given in Fig. 3.2. This integration strategy appears as a very interesting solution in many aspects. First of all, if less than four satellites are tracked the DR errors can still be estimated, so this approach is highly suitable for urban environments. Moreover, fault detection techniques can be applied to each GPS raw measurement to prevent positioning errors introduced by outliers [START_REF] Kubrak | Hybridisation of a GPS receiver with low-cost sensors for personal positioning in urban environments[END_REF]. However, it presents the disadvantage of introducing non linear measurement equations, so that its implementation is more complex and expensive than the loosely coupled strategy.

Chosen integration

Before stating the integration strategy adopted in this work, some remarks must be done on the tight and loose integration performances in urban scenarios. As already explained, tight coupling is generally considered as the best urban navigation solution because of its tolerance to limited visibility over the GPS constellation. Even in the presence of less than 4 GPS satellites the INS can still be corrected. However, a main problem arises when the GPS measurement to provide inertial corrections are subjected to important errors. As corrections are done at a higher level in the tight integration approach, a non detected outlier will have a stronger impact on the position solution than when using a loose integration strategy where the already filtered GPS position solution is compared to the DR-based position. Of course, the more accurate the DR sensors, the lower the impact of the outlier.

Nevertheless, this presents an import trade off to be considered when performing urban navigation with relative low quality sensors. On the other hand, as shown in [START_REF] Petovello | Real-Time Integration of a Tactical-Grade IMU and GPS for High-Accuracy Positioning and Navigation[END_REF], no significant differences are found between the loose and tight integration in open sky scenarios.

The main idea of this chapter is to introduce a suitable integration strategy that will be further exploited in the next chapters to develop fault detection techniques for GPS signals in urban scenarios.

Hence, tight integration is chosen as the coupling approach to be used herein. As already explained, one of the main drawbacks of the tight integration approach is the use of non linear measurement equations for the pseudorange or Doppler measurements. An Unscented Kalman Filter (UKF) will be implemented to test the impact of such non linearities. Results will be compared to those obtained when using a standard Extended Kalman filter (EKF).

Filtering framework

The estimation principle for dynamic systems consist of recursively evaluating the evolution of the system variables X from some available noisy observations Y. In particular, we are going to work with 66CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA markovian systems that can be modeled by the following state-space equations:

X t = f t (X t-1 ) + v t , (3.1)
Y t = h t (X t ) + w t . (3.2)
Eq. (3.1) describes the a priori dynamics followed by the unknown parameters X t and eq. 3.2 is the measurement model that associates the parameters with the measurements Y t , where

• X t ∈ R nx is the state vector containing the parameters to be estimated,

• f t is the non linear state transition function,

• Y t ∈ R ny is the measurement vector (also called observation vector),

• h t is the non linear measurement function,

• v t ∈ R nx and w t ∈ R ny are additive white noises.

The filtering problem is resolved by calculating an estimator that maximises a given criterion that it is usually application dependent. The filtering techniques presented in the following applied a minimum mean square error (MMSE) criterion. The MMSE estimator ( Xt ) MMSE is expressed as:

( Xt ) MMSE = E [X t |Y 1:t ] = Xt X t p(X t |Y 1:t )dX t , (3.3)
where Y 1:t = (Y 1 , . . . , Y t ) and p(X t |Y 1:t ) is the conditional probability density (also called a posteriori probability density) of the state vector considering all the received measurements. The dynamic nature of the above described system implies a growing number of measurements Y 1:t with time. To avoid having an increasing system complexity, a recursive approach is used to estimate the probability density at every time instant t. This approach is separated in two different stages:

1. Prediction

The state vector is propagated at time instant t conditional to the previous received measurements and according to an a priori knowledge on the state evolution. From the Chapman-Kolmogorov equation we obtain:

p(X t |Y 1:t-1 ) = X t-1 p(X t |X t-1 )p(X t-1 |Y 1:t-1 )dX t-1 , (3.4)
where p(X t |X t-1 ) is the state transition prior.

Correction

The predicted state is updated by the received measurement according to Bayes rule:

p(X t |Y 1:t ) = p(Y t |X t )p(X t |Y 1:t-1 ) p(Y t |Y 1:t-1 ) , (3.5)
where p(Y t |X t ) is the observation likelihood density, and the normalizing factor p(Y t |Y 1:t-1 ) is obtained from

p(Y t |Y 1:t-1 ) = Xt p(Y t |X t )p(X t |Y 1:t-1 )dX t . (3.6)
The application of these two stages provides the optimal filtered solution for a non linear system.

In practice, eq. (3.4) and (3.6) involve the calculation of multidimensional integrals that can be difficult to evaluate. However, for the specific case of linear systems with Gaussian random variables the optimal solution can be computed, and it is given by the well known Kalman filter proposed in [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] for the discrete case and in [START_REF] Kalman | New results in linear filtering and prediction theory[END_REF] for the continuous case. For most real world models (such as the navigation model addressed in this thesis) the multi-dimensional integrals are intractable and approximated solutions must be used. The Extended Kalman Filter (EKF) and Unscented Kalman filter (UKF) will be discussed with this purpose.

The optimal Kalman Filter (KF)-Linear systems

Kalman [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF] derived a recursive form to update the conditional mean of the state random variables Xt = E[X t |Y 1:t ] (i.e. to compute MMSE estimator given in eq. (3.3)) and its covariance P t , given by the following expressions:

Xt = Xt|t-1 + K t Y t -Ŷt|t-1 , (3.7)
P t = P t|t-1 -K t P t|t-1 K T t , (3.8)
where Xt|t-1 is the predicted state vector, Ŷt|t-1 is the predicted measurement vector, and K t is the optimal gain matrix. The gain matrix is one of the crucial points in the Kalman filter strategy and it can be understood as a weighting factor indicating how much of the new information contained in 68CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA the measurements should be accepted by the system. Generally speaking, it weights the information from the measurements against the current knowledge of the states.

While the recursion presented in (3.7)-(3.8) is linear, it must be noted that no linearity of the model was assumed. At the same time, no assumptions are made on the noise distribution or its Gaussian behavior. In fact, Kalman stated that consistent minimum variance estimates of the system random variables and hence the posterior state distribution can be calculated by maintaining only the first and second moments. In this way, general expressions for the predicted state variables, measurements and Kalman gain are given by:

Xt|t-1 = E [f t (X t-1 ) + v] , (3.9
)

Ŷt|t-1 = E [h t (X t ) + w] , (3.10)
K t = E (X t -Xt|t-1 )(Y t -Ŷt|t-1 ) T E (Y t -Ŷt|t-1 )(Y t -Ŷt|t-1 ) T -1 (3.11) = P XY P -1 Ŷ , (3.12)
where P XY is the expected cross-correlation matrix (covariance matrix) of the state prediction error and the observation predicted error, and P Ŷ is the expected auto-correlation matrix of the observation predicted error. Indeed, for the optimal solution to be computed (i.e to obtain exact "analytical" results for (3.9)-(3.12)) the filter must be applied to linear Gaussian systems. Under this condition, the Bayesian recursion in 3.5 can be exactly evaluated, and the computed mean and covariance parameters (3.7)-(3.8) fully characterized the resulting Gaussian noise. The dynamic state-space model for a linear Gaussian system is written as:

X t = F t X t-1 + v t , (3.13)
Y t = H t X t + w t , (3.14)
where F t and H t are the linear transition and measurement functions. The noise terms v t and w t follow zero-mean Gaussian distribution with v t ∼ N (0, Q t ) and w t ∼ N (0, R t ), where Q t and R t are the noise covariance matrices. The prediction and correction steps are thus described by the following equations:

Prediction

State propagation

Xt|t-1 = F t Xt-1 (3.15)

State covariance matrix propagation

P t|t-1 = F t P t-1 F T t + Q t (3.16)
Predicted measurement from the state vector variables

Ŷt|t-1 = H t Xt|t-1 (3.17)

Correction

Kalman's gain

K t = P t|t-1 H T t H t P t|t-1 H T t + R t -1 (3.18)
State vector update

Xt = Xt|t-1 + K t Y t -Ŷt|t-1 (3.19)
State covariance matrix update

P t = (I nx -K t H t ) P t|t-1 (3.20)
where the error between the true observation and the predicted observation Y t -Ŷt|t-1 is usually called the innovation, whose covariance matrix is given by P Ŷ = H t P t|t-1 H T t + R t (P Ŷ will be also referred hereafter as S t ). More details about the Kalman filter can be found in [START_REF] Anderson | Optimal filtering[END_REF].

Non-linear systems

As it was already explained, non linear systems prevent the Bayesian recursion to be analytically resolved. Therefore, suboptimal strategies must be implemented to address this problem. In particular two categories summarize the different types of approximations that are generally applied:

• the non linear function(s) are approximated by their linearized first order Taylor expansion, 70CHAPTER

HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

• the approximation is done on the a posteriori distribution p(X t |Y t ) (i.e. just a limited number of moments is used to describe the distribution).

In particular, the Extended Kalman Filter is based on the first proposal, while the Unscented Kalman Filter addresses the second proposal.

The Extended Kalman Filter (EKF)

The extended Kalman filter, detailed in [START_REF] Anderson | Optimal filtering[END_REF] and [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF], is an "extension" of the Kalman filter to non linear dynamic state-space models. It linearizes the system around the current estimate using a first order truncation of the multi-dimensional Taylor series expansion. If the state variables do not exhibit strong dynamics withing the sampling period, the linearized transformations can be considered reliable (i.e. no important information is lost). The a posteriori distributions are still approximated by Gaussian posterior density functions (pdfs). The first and second order moments (i.e. mean and variance) are propagated using the optimal KF equations. The linearized state model is constructed as follows

X t ≈ F t X t-1 -X t-1 + f t X t-1|t-1 , (3.21)
with

F t = δf t δX t-1 Xt-1|t-1 . (3.22)
Similarly, the linearized measurement model is written as

Y t ≈ H t X t-1 -X t|t-1 + h t X t|t-1 , (3.23)
with

H t = δh t δX t Xt|t-1 . (3.24)
F t and H t are usually referred as the Jacobian matrices of f t and h t respectively. The full recursive algorithm is given in Table 3.1.

• Initialization X0 = E [X 0 ] , P 0 = E X 0 -X0 X 0 -X0 T • For t = 1, . . . , end:
1. Prediction

• Compute the state linearized transition matrix

F t = δft δX t-1 δX t-1|t-1
• Compute the predicted state mean and covariance

X t|t-1 = f t X t-1 P t|t-1 = F t P t-1 F T t + Q t • Compute the linearized measurement matrix H t = δht δXt δX t|t-1

• Compute the predicted observation mean, innovation and innovation covariance

Ŷt|t-1 = h t Xt|t-1

I t = Y t -Ŷt|t-1 S t = H t P t|t-1 H T t + R t

Correction

• Compute Kalman's gain

K t = P t|t-1 H T t (S t ) -1

EKF limitations

The "first order Taylor series linearization" used by the EKF introduces a number of serious limitations to its application domain. A brief discussion of its most important flaws is given in the following. An illustration of the distortions introduced by the linearization process is given in Fig. 3.3.

1. The linearized transformations are reliable if the error propagation can be correctly approximated by a linear function during the sampling period. Otherwise, it can lead to important problems in the stability of the filter [START_REF] Giremus | Apports des techniques de filtrage non linéaire pour la navigation avec les systèmes de navigation inertiels et le GPS[END_REF].

2. The computation of the linearized or jacobian matrices cannot always be possible. Systems can be subjected to abrupt changes, singularities or discontinuities that would unable the jacobian computation. At the same time, depending on the complexity of the system, the calculation of the jacobians can entail an unfordable large computational cost.

3. As explained in [vdM04], the first order Taylor series linearization ignores the fact that the state variables are random variables themselves by failing to account for the "probabilistic spread" of these variables. It is important to understand that though the state covariance matrix P t plays an important role in the validity of the EKF "first order linearization", it is completely ignored during the linearization process. In other words, when the predicted measurement mean Ŷt|t-1 is computed just as a function of the predicted state mean Ŷt|t-1 = h t Xt|t-1 , no higher orders of the Taylor series expansion accounting for the actual variable spread (X t -Xt|t-1 )

are considered. This linearization around a single point can have serious consequences on the accuracy and consistency of the filter.

The unscented Kalman filter [START_REF] Julier | Unscented Filtering and Nonlinear Estimation[END_REF] was developed to address the deficiencies of the EKF linearization by providing a derivativeless strategy using several propagated points that enable a correct characterization of the state a posteriori mean and covariance up to a second order.

Nonlinear transformation

Linearized transformation PDF PDF

Real

Figure 3.3: Linearized transformation (EKF).

The Unscented Kalman Filter (UKF)

The unscented Kalman filter introduced by Julier and Uhlmann in [JU96; JU04] is aimed at overcoming the flaws presented by the EKF due to the need of linearizing the system equations. It is a minimum mean square (MMSE) estimator based on the optimal Kalman filter framework. In the UKF the state distribution is still represented by a Gaussian random variable, but it is specified by a minimal set of point (sigma points) which are deterministically chosen. These points completely capture the true mean and covariance of the state vector variables. So when propagated through the true non-linear system, they provide a posterior mean and covariance accurately to the second order of any non linearity. In comparison, the EKF only calculates the posterior mean and covariance accurately to the first order by truncating all the higher orders. The criterion for choosing the adequate set of sigma points is based on the unscented transformation [START_REF] Julier | A general method for approximating non-linear transformations of probability distributions[END_REF]. The idea is to overcome the EKF linearization drawbacks by not incuring in the implementation nor computational cost of higher order filtering schemes. The principles for this transformation are presented in the next section.

Unscented transform (UT)

The UT is built on the principle that it is easier to approximate a probability distribution than an arbitrary non linear function. The method allows to calculate the statistics of a random variable that undergoes a non linear transformation by using a set of weighted sigma points [START_REF] Julier | A general method for approximating non-linear transformations of probability distributions[END_REF]. The following steps summarize the UT approach:

1. A set of sigma points χ i and their associated weights W i is computed so at least the mean 74CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA X and covariance P of the prior random variables are totally captured. The weights can be positive or negative, but to provide an unbiased estimate they must obey the condition:

2nx i=0 W i = 1. (3.25)
where n x is the dimension of vector X.

2. The sigma points are propagated through the true non linear function h(•) to generate a set of transformed sigma points according to:

Υ i = h (χ i ) . (3.26)
3. The posterior statistics (i.e. mean, covariance and cross-covariance) are computed by using the propagated sigma points and their respective weights as follows: Note that Ŷ , P Ŷ and P XY are all second order approximations. The UT principle is illustrated in Fig. 3.4. A set of weighted sigma points that satisfies the above mentioned conditions consists of a symmetric set given by: point locations such that the first and second order terms match exactly. Details on the sigma-point computation can be found in [vdM04]. The scaling factor λ is aimed at providing an additional degree of freedom to control some aspects of the higher moments of the distribution of the sigma points. It was discussed in [START_REF] Julier | A general method for approximating non-linear transformations of probability distributions[END_REF] that λ = 3 -n x is optimal for a single-state Gaussian distribution, and for multi-dimensional systems, if 0 < n x + λ < 3, then the absolute error in the predicted mean is smaller than with linearization.

Ŷ = 2nx i=0 W i Υ i , (3.27)
P Ŷ = 2nx i=0 W i (Υ i -Ŷ )(Υ i -Ŷ ) T , (3.28)
P XY = 2nx i=0 W i (χ i -X)(Υ i -Ŷ ) T , ( 3
χ 0 = X, W 0 = k (n x + k) (3.31) χ i = X + (n + k)P i , W i = 1 2(n x + k) (3.32) χ i+nx = X - (n x + k)P i , W i+nx = 1 2(n x + k) ( 3
It is important to understand that though the UT appears to share common principles with the particle filter, there are several fundamental differences. First, the sigma points are not randomly but deterministically chosen. Second, the associated weights W i can take values that lie outside the [0, 1]

interval, which will be considered inconsistent with the particle filter sampling approach. Moreover, the computational cost of the UT algorithm is the same order of magnitude as the EKF.

The sigma-point selection used in the UT has the property that as the dimension of the state-space (n x ) increases, the radius of the sphere containing the sigma points increases as well. Though the mean and covariance of the prior distribution are still correctly captured, as the χ i,i =0 sigma points get spread further away from the mean of the distribution, the possibility of sampling undesired non-local effects is important. If the nonlinearities in question are important, this can lead to severe difficulties [vdM04]. The scaled unscented transformation proposed in [START_REF] Julier | The scaled unscented transformation[END_REF] helps to address this problem by scaling the sigma points to be used.

The scaled unscented transformation (SUT)

An interesting example is given in [START_REF] Shin | Estimation techniques for low-cost inertial navigation[END_REF] to understand the necessity for a scaled unscented tranformation (SUT) when working with low-cost INSs. For instance, imagine that the heading angle error from a highly inaccurate gyroscope must be estimated. Consider a high dimension state vector with a spreading factor for the sigma points (n x + k) ≈ 4, and a heading angle uncertainty of 30 • . The heading angle will thus be sampled in the range of ±120 • enabling non-local effects to be introduced in the transformation. Furthermore, attitude angles do not belong to a vector space in the sense they have a period of 360 • and they thus repeat themselves. Hence, for large heading uncertainties, some angles can be sampled twice. To avoid this effect, the SUT [START_REF] Julier | The scaled unscented transformation[END_REF] replaces the original set of sigma points with a transformed set given by

χ i = χ 0 + α(χ i -χ 0 ) i = 0, . . . , 2n x , (3.34)
where α is a positive scaling parameter that can be made arbitrary small (0 < α < 1) to minimize the sampling of non local effects. The new set of sigma point weights is then computed as:

W i =    W 0 /α 2 + (1 -1/α 2 ) i = 0 W i /α 2 i = 1, . . . , 2n x (3.35)
where indeed, a further parameter β will be introduced to allow for the minimization of higher order errors if prior knowledge (e.g. kurtosis) of the distribution of X is available. The zeroth sigma point χ 0 will be therefore affected with two types of weights: one used for the mean calculation in (3.27) and the other one for the covariance calculations in [(3.28),(3.29)] which is defined as a function of β.

The final set of sigma points is thus described as:

χ 0 = X,    W (m) 0 = λ (nx+λ) W (c) 0 = λ (nx+λ) + (1 -α 2 + β) (3.36) χ i = X + (n x + λ)P i , W (m) i = W (c) i = 1 2(n x + λ) (3.37) χ i+nx = X - (n x + λ)P i , W (m) i+nx = W (c) i+nx = 1 2(n x + λ) (3.38)
where λ = α 2 (L + k) -L and β is a non negative weighting term (β = 2 is optimal for Gaussian distributions [START_REF] Julier | Unscented Filtering and Nonlinear Estimation[END_REF]). W covariance calculation.

Filter implementation

The UKF is a direct application of the unscented transformation to the recursive Kalman filter.

However, to directly apply the Kalman filter principles to the sigma points obtained in (3.36)-(3.38), the state vector is generally augmented with the process and noise terms such as:

X a t =      X x t X v t X w t      =      X t v t w t      (3.39)
where X a t is the new augmented state vector of dimension n a x = 2n x + n y and n y is dimension of the measurement vector. The sub-vectors X x t , X v t and X w t refer to the nominal EKF state vector, state noise and measurement noise terms, respectively. Similarly, the new augmented state covariance matrix is given by

P a t =      P t 0 nx×nx 0 nx×ny 0 nx×nx Q t 0 nx×ny 0 nx×nx 0 nx×nx R t      (3.40)
where P t is the same state covariance matrix used for the EKF, and Q t and R t are the state and measurement noise covariance matrices. By including the noise terms in the augmented system, the noise statistics are captured with the same accuracy as the state variables. The UKF algorithm is described in table 3.2, where χ a = (χ x , χ b , χ w ) T and Xt and P t denote the final estimated mean and covariance matrix of the nominal state vector X t .

In particular, when the state and measurement noises are additive (zero mean) noises, the complexity of the previous algorithm can be reduced. The noise terms need not to be included in the state vector, so both the quantity and dimension of the sigma points are decreased [START_REF] Wan | The Unscented Kalman Filter[END_REF][START_REF] Van Der Merwe | Sigma-Point Kalman Filters for Probabilistic Inference in Dynamic State-Space Model[END_REF].

Though the main structure of the UKF algorithm given in table 3.2 remains unchanged, some expressions need to be modified to account for the noise contribution. Table 3.3 summarizes the UKF algorithm in the presence of additive zero mean noise.

• Initialization X0 = E [X 0 ] , P 0 = E X 0 -X0 X 0 -X0 T Xa 0 = E[X a 0 ] = [ X0 , 0 1×nx , 0 1×ny ] T , P a 0 =     P 0 0 nx×nx 0 nx×ny 0 nx×nx Q 0 0 nx×ny 0 nx×nx 0 nx×nx R 0     ,
• For t = 1, . . . , end:

• Calculate sigma points and weights: • Time-update sigma points:

χ a t-1 = Xa t-1 , Xa t-1 + (n a x + λ)P a t-1 , Xa t-1 -(n a x + λ)P a t-1 W a t-1 = W 0 , W i , W i+n a x ,
χ x t|t-1 = f(χ x t|t-1 , χ v t|t-1 )
• Calculate mean and covariance of the predicted sigma points:

Xt|t-1 = 2n a x i=0 W a,(m) i χ x i,t|t-1 P t|t-1 = 2n a x i=0 W a,(c) i χ x i,t|t-1 -Xt|t-1 χ x i,t|t-1 -Xt|t-1 T
• Apply the observation model to the predicted sigma points:

Υ i = h χ x t|t-1 , χ w t|t-1
• Calculate mean, covariance and cross-covariance associated to the predicted observations:

Ŷt = 2n a x i=0 W a,(m) i Υ i,t|t-1 P Ŷt = 2n a x i=0 W a,(c) i Υ i,t|t-1 -Ŷt Υ i,t|t-1 -Ŷt T P XtYt = 2n a x i=0 W a,(c) i χ x i,t|t-1 -Xt|t-1 Υ i,t|t-1 -Ŷt T
• Update the sate vector using standard KF equations:

K t = P XtYt P -1 Ŷt I t = Y t -Ŷt Xt = Xt|t-1 + K t I t P t = P t|t-1 -K t P Ŷt K T t Table 3.2: UKF algorithm. • Initialization X0 = E [X 0 ] , P 0 = E X 0 -X0 X 0 -X0 T • For t = 1, . . . , end:
• Calculate sigma points and weights:

χ t-1 = Xt-1 , Xt-1 + (n x + λ)P t-1 , Xt-1 -(n x + λ)P t-1 W t-1 = [W 0 , W i , W i+nx ] , from eqs. (3.36)-(3.38) with i = 1, ..., n x .
• Time-update sigma points:

χt|t-1 = f(χ t|t-1 )
• Calculate mean and covariance of the predicted sigma points:

Xt|t-1 = 2nx i=0 W i χi,t|t-1 P t|t-1 = 2nx i=0 W i χi,t|t-1 -Xt|t-1 χi,t|t-1 -Xt|t-1 T + Q t
• New set of redistributed sigma points to account for state noise:

χ t|t-1 = χt|t-1 , χ0,t|t-1 + (n x + λ)Q t , χ0,t|t-1 -(n x + λ)Q t
• Apply the observation model to the predicted sigma points:

Υ i = h χ t|t-1
• Calculate mean, covariance and cross-covariance associated to the predicted observations:

Ŷt = 2nx i=0 W (m) i Υ i,t|t-1 P Ŷt = 2nx i=0 W (c) i Υ i,t|t-1 -Ŷt Υ i,t|t-1 -Ŷt T + R t P XtYt = 2nx i=0 W (c) i χ i,t|t-1 -Xt|t-1 Υ i,t|t-1 -Ŷt T
Expressions for K t , I t , Xt and P Xt are the same as for 3.2.

Cramer Rao bound

All the filtering approaches presented so far provide suboptimal approximated solutions to the navigation problem. To analyze the quality of the estimations and to quantify the errors introduced by approximations, it is of great practical interest to compare the computed quadratic errors to a lower bound. Within this context, the Cramer-Rao bound is normally used because it provides a lower bound on the estimator mean square error [START_REF] Jazwinski | Stochastic processes and filtering theory[END_REF]. The Cramer Rao bound can be viewed as a reference to which the state MSEs of suboptimal algorithms can be compared. It is a mathematical statistical tool that can be applied to any kind of linear or non linear statistical system with Gaussian or non Gaussian noise. In particular, the additive Gaussian noise will be studied in this thesis.

The bound is strictly related to the information introduced in the system by the observations. It is obtained as the inverse of the Fisher information matrix. For a linear Gaussian case, the bound and the KF covariance matrix are the same.

This thesis will focus on the posterior Cramer-Rao bound (PCRB) which is often referred to as the Bayesian version of the Cramer-Rao bound [START_REF] Van Trees | Detection, estimation and modulation theory[END_REF]. Usually, the Cramer Rao bound is defined in the parametric framework and it is based on the likelihood of the estimation problem. It is applied to unknown but fixed parameters. On the other hand, the PCRB is aimed at random parameters which are represented by their prior distribution. Here, the joint density of parameters and observations plays the role of the likelihood in the parametric bound. Another difference, resides in the fact that the PCRB can be applied to biased estimators where the a priori conditions for the unknown parameter distributions are weak. Recently, in [START_REF] Tichavsky | Posterior Cramer-Rao bounds for-discrete time nonlinear filtering[END_REF] and [START_REF] Bergman | Posterior Cramer-Rao Bounds for Sequential Estimation[END_REF], extensions to the use of the PCRB to nonlinear discrete systems were proposed.

Considering a non-linear system given by eqs. (3.1) and (3.2), the Fisher information matrix J t at time instant t is given by Therefore, for any estimator X(Y 1:t ), the PCRB verifies

J t = E p(Xt,Y 1:t ) - ∂ 2 log p(X t ,
P t = E ( X(Y 1:t ) -X t )( X(Y 1:t ) -X t ) T ≥ J -1 t = P CRB t (3.43)
where the matrix comparison given by the logic operator ≥ means that the difference between the expected mean square error correlation matrix and the PCRB matrix is positive semi-definite.

All the filtering tools presented so far will be applied to different configurations of hybrid GPS/DR land navigation systems. The performances of the EKF and UKF will be compared to the optimal CRB. The validity of the filtering assumptions done for both the EKF and UKF will be analyzed.

GPS/Dead-reckoning-sensors hybrid systems

This section will present two multi-aided land navigation systems combining either GPS and WSS measurements (GPS/WSS), or further including the INS outputs (GPS/INS/WSS).

GPS/WSS integration

The parameters to be estimated by the filter include the WSS and GPS errors. Indeed, GPS measurements will aid to update WSS errors. The error model for the WSS measurements was already introduced at the end of chapter 2. Just the two main equations describing the measured velocity and yaw rate are recalled herein,

V t = ω rl t R rl + ω rr t R rr 2 + n V t . ( 3
δR rr t = n r t , δR rl t = n l t ,
where n r t and n l t are white noise sequences.

82CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

State model

The wheel radius is going to be considered the main source of error affecting the WSS measurements and it will be thus included in the filter state vector. Similarly, GPS errors must also be modeled and predicted. In chapter 1, several error sources affecting the GPS measurement were described. In particular, errors induced by satellite-based inaccuracies (i.e. satellite clock error or ephemeris errors)

and atmospheric propagation, such as the ionospheric and tropospheric error, are considered to be correctly compensated by EGNOS. Among the remaining errors, multipath or interferences are not usually included in the state vector because they are application dependent and therefore difficult to be modeled. Therefore, just the receiver clock bias and drift are estimated. A second order model adequately describes the error behavior:

  ḃt ḋt   =   0 1 0 0   F GP S t   b t d t   +   σ b 0 0 σ d   BGP S t w t , (3.46)
where b t denotes the GPS receiver clock offset in meters, d t is the derivative of b t and w t is a Gaussian noise vector w ∼ N (0, I 2 ). The standard deviations of the noise distributions (σ b , σ d ) depend on the type of local oscillator used by the receiver [START_REF] Kubrak | Hybridisation of a GPS receiver with low-cost sensors for personal positioning in urban environments[END_REF]. The resulting time-continuous state vector for the GPS/WSS strategy is composed of the errors in the vehicle position, velocity, angle of direction (heading), radii and receiver clock parameters:

Ẋt = δ λt , δ φt , δ
F GP S/W SS t =   F W SS t 0 0 F GP S t   , BGP S/W SS t =   BW SS t 0 0 BGP S t   . (3.48)

Measurement model

The GPS/WSS hybrid strategy uses the GPS measurements to correct the WSS navigation errors. In particular the pseudorange measurements will be exploited in this chapter. The observation vector Y t will be then composed by the pseudoranges corresponding to the n y visible satellites:

Y t = ρ 1 , . . . , ρ ny T (3.49)
where ρ i denotes the ith pseudorange. The pseudorange expression was already given in chapter 1

and it is recalled in eq. (3.50):

ρ i = (x -x s i ) 2 + (y -y s i ) 2 + (z -z s i ) 2 + b t + n ρ i (3.50)
where i = 1, . . . , n y , (x, y, z) is the vehicle position and (x s i , y s i , z s i ) is the ith satellite position both in the in the ECEF frame. The measurement equation is thus obtained by relating the state vector variables in (3.47) to the pseudorange expression in (3.50). A frame transformation is therefore necessary to convert the geodetic coordinates (λ W SS t + δλ t , ϕ W SS t + δϕ t , h t ) into cartesian coordinates (x, y, z):

g(δλ t , δϕ t )          x t = (N + h t ) cos(λ W SS t + δλ t ) cos(ϕ W SS t + δϕ t ), y t = (N + h t ) cos(λ W SS t + δλ t ) sin(ϕ W SS t + δϕ t ), z t = [N (1 -e 2 ) + h t ] sin(λ W SS t + δλ t ) (3.51) where N = a √
1-e 2 sin(λ) 2 is the length of the normal to the ellipsoid, from the surface of the ellipsoid to its intersection with the ECEF z-axis frame, for a = 6378137m (ellipsoid semimajor axis length) and e = 0.0818 (eccentricity of the ellipsoid). For the assumed 2D model, height h t is considered known (i.e. it will be given by the high quality navigation reference system described in section 3.8.2). The ECEF position coordinates can be then used to calculate the n y received pseudoranges according to the following measurement model [START_REF] Giremus | Apports des techniques de filtrage non linéaire pour la navigation avec les systèmes de navigation inertiels et le GPS[END_REF]:

Y t = h t (X t ) + n ρ , n ρ ∼ N (0, R t ) (3.52) h t (X t ) = q t • g(δλ t , δϕ t ) + b t = q t (g(δλ t , δϕ t )) + b t (3.53) 84CHAPTER 3.

HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

where a • b denotes the function composition a composed with b, and q t is the function for the range measurements: Assuming the pseudoranges are independent, the measurement noise covariance matrix R t is expressed as

q t (x t , y t , z t ) =              x t -
R t =         σ 2 t,1 0 • • • 0 0 . . . . . . . . . . . . 0 0 • • • 0 σ 2 t,ny         . (3.55)
where the pseudorange variances are a function of the carrier to noise power ratio C/N 0 measured by the receiver and the EGNOS corrections as explained in chapter 1.

GPS/INS/WSS integration

The commonly used state-space model for the GPS/INS hybrid system will be first introduced. Based on its nominal dynamic state-space model the system will be further augmented to include navigation data provided by the WSS outputs. Indeed, the "augmented" state vector will include the wheel radius errors while the observation model will include the vehicle along track velocity and yaw angle computed from the WSS outputs. The idea is that the WSSs will enhance the system performances by providing additional measurements for the estimation of the INS errors. Indeed, in the presence of GPS signals, the WSS contribution will not be significant [START_REF] Gao | Development of precise GPS/INS/Wheel Speed Sensor/yaw rate sensor integrated vehicular positioning system[END_REF]. However, during GPS partial or total outages, WSS measurements will help to control the growing INS errors.

GPS/INS

The GPS/INS hybrid system is one of the most well-known integration techniques. Several references can be found in the literature concerning the performances of GPS/INS systems applied to different scenarios. The common GPS/INS coupling techniques use the GPS pseudorange measurements to correct the INS errors in order to obtain reliable position estimation. For instance, [START_REF] Petovello | Real-Time Integration of a Tactical-Grade IMU and GPS for High-Accuracy Positioning and Navigation[END_REF] studied how the integration of a high grade tactical IMU significantly improves both the solution accuracy and the system reliability. By reliability, we mean the capacity to detect outliers when the received measurement is compared to a given bound or protection level. The inclusion of a good quality IMU to the navigation system helps the uncertainty in the position estimation to be reduced and therefore fault/outlier detection performances are enhanced.

The variables to be estimated for the GPS/INS system generally include the different sources of errors affecting the GPS and INS measurements. The time-continuous state vector is thus constructed as follows,

Ẋt = δ ṗt , δ vt , δ ψt , ḃa,t , ḃg,t , ḃt , ḋt T GPS/INS ∈ R 10 , (3.56)
where (b t , d t ) describe the GPS clock errors affecting the pseudorange measurements, and the subvector (δ ṗt , δ vt , δ ψt , ḃa,t , ḃg,t ) contains the INS errors given by

• δ ṗt = δ λt , δ φt denotes the INS position errors in the navigation frame,

• δ v = (δ vn,t , δ v e,t ) denotes the INS velocity errors in the navigation frame,

• δ φt is the heading (or yaw) angle error,

• ḃa,t is the two dimension accelerometer bias and ḃg,t is the gyro bias associated to the heading angle.

The

F GP S/IN S/W SS t =      F GP S/IN S t 0 0   0 0 0 0        , BGP S/IN S/W SS t =      BGP S/IN S t 0 0   σ rl 0 0 σ rr        , (3.59)
where σ rl and σ rr are the noise standard deviations for the left and right wheel radius errors.

Measurement model

The proposed GPS/INS/WSS integration uses not only the pseudorange measurements, but also the along track velocity and yaw rate computed from the WSS as measurements to correct the INS errors. In particular, non-holonomic constraints are going to be considered for the INS velocity vector in the mobile frame as in [START_REF] Gao | Development of precise GPS/INS/Wheel Speed Sensor/yaw rate sensor integrated vehicular positioning system[END_REF]. The idea is to derive the observations from constraint equations reflecting the behavior of a land vehicle so as to have a maximum number of measurements available during GPS outages. Assume that the side slip effects can be neglected as in [START_REF] Enquist | Aspects of high precision estimation of vehicle dynamics[END_REF], and there is no misalignment between the body frame and the vehicle. Therefore, it is valid to say that the cross-track velocity should be zero v y = 0 (i.e. using a non-holomonic constraint [START_REF] Gao | Development of precise GPS/INS/Wheel Speed Sensor/yaw rate sensor integrated vehicular positioning system[END_REF]). These additional virtual measurements will mainly help to keep track of INS errors under degraded GPS visibility.

Indeed, for civil land vehicle navigation, the cross-track accelerations are usually so weak that they are almost masked by the low cost sensor noise. Within this context, better performances can be generally obtained when setting the velocity constraint v y = 0 [God06; Shi05] (a noise component is also associated to this constraint as later presented in (3.76)). However, a frame transformation has to be used to compare the velocities given by the WSS and INS. Here, INS velocities are transformed into the vehicle frame (an alternative would be to express WSS velocities in the navigation frame).

The WSS velocity measurement equation is given by • n vm t ∈ R 2 is a white noise vector n vm t ∼ N (0, R v t )

  V t 0   WSS -R T m2n (v n t ) INS = h v (
• R T m2n denotes the navigation to mobile transformation matrix as defined in chapter 2 (from the equivalence R n2m = R T m2n ),

• (v n t ) INS is the vehicle velocity in the navigation frame obtained from the INS.

The expression for h v (.) is given by

h v (δR rr t , δR rl t , δψ t , δv n t ) = δ(R T m2n (v n t ) INS ) -   δV t 0   WSS (3.61)
where the error propagation was done according to a + δa = b + δb ⇒ a -b = δb -δa. The expression for δ(R T m2n (v n t ) INS ) is obtained by propagating the errors present in R T m2n and v n t as

δ(R T m2n (v n t ) INS ) = δ(R T m2n )(v n t ) INS + R T m2n δ(v n t ) INS = R T m2n V n INS,t δψ t + R T m2n δ(v n t ) INS (3.62)
where the term R T m2n V n INS,t δψ t is derived from the small angle transformation assumption (as it was implicitly done when presenting the INS error model in chapter 2) [START_REF] Farrell | The Global Positioning System and Inertial Navigation[END_REF]. The principles for this assumption are explained in the following. Consider R * m2n the real 2D transformation matrix from frame m to frame n and R m2n the analogous transformation matrix but computed with the heading angle given by the INS. R m2n will differ from R * m2n in a additional rotation caused by the error in the INS heading angle δψ t . The two transformation matrices can be thus related as follows:

R * m2n =   cos(δψ t ) -sin(δψ t ) sin(δψ t ) cos(δψ t )   R m2n . (3.63)
Consider now a very small angle error so that cos(δψ t ) and sin(δψ t ) can be correctly approximated by cos(δψ t ) ≈ 1 and sin(δψ t ) ≈ δψ. In that case, (3.63) can be rewritten as

R * m2n =   1 -δψ t δψ t 1   R m2n (3.64) R * m2n =        I 2 +   0 -δψ t δψ t 0   ξ        R m2n (3.65) = R m2n + R m2n ξ = R m2n + δR m2n
In particular, as we are interested in the navigation n to mobile m transformation R n2m , it is useful to introduce the following relation [FB99, p.200]:

R * n2m = (R * m2n ) T = R n2m   I 2 -   0 -δψ t δψ t 0     (3.66)
where this relation is accurate to first order and can be easily verified by considering (I -ξ) -1 = (I + ξ) to first order and ξ T = -ξ. Therefore, the term associated to the error introduced by the transformation matrix in (3.62) is expressed as:

δ(R T m2n )(v n t ) INS = R T m2n (-ξ) T (v n t ) INS (3.67)
where for convenience the right hand side of (3.67) is explicitly expressed as a linear function of the state vector variable δψ t as follows

δ(R T m2n )(v n t ) INS = R T m2n   0 -δψ t δψ t 0     v n v e   (v n t ) INS = R T m2n   -v e v n   V n INS,t δψ t (3.68)
The last term in the right hand of eq. (3.61) represents the WSS velocity error in terms of the radii errors. This term can be easily deduced from eq.(3.44) by considering that the actual wheel radius differs from the nominal radius R in δR. The expression for δV t , which was already given in chapter 2 and used to construct the F W SS t matrix used in (3.48), is recalled here:

δV t = ω rl t δR rl t + ω rr t δR rr t 2 . (3.69)
The measurement model for the WSS yaw angle is straightforward since no frame transformation must be done. Similarly to (3.60), the measurement model is given by the difference between the two DR navigation systems as:

(ψ t ) WSS -(ψ t ) INS = h ψ (
Y a t =               ρ 1 . . . ρ ny V t -(v x t ) INS -(v y t ) INS (ψ t ) WSS -(ψ t ) INS               T (3.72)
where (v x t , v y t ) INS are the INS velocities in the mobile frame. Hence, the augmented measurement model is described as:

Y a t = h a (X t ) GPS/INS/WSS + n a , (3.73)
h a (X t ) GPS/INS/WSS =      h(•) h v (•) h ψ (•)      , (3.74)
where

• h a (•) represents the augmented measurement function formed by the concatenation of h(•), h v (•), and h ψ (•),

• h(•) is the measurement function for the pseudorange measurements described in (3.52) where

(λ W SS t , ϕ W SS t ) are replaced by (λ IN S t , ϕ IN S t
),

• h v (•) is the measurement function based on the WSS velocity measurement described in (3.61),

• h ψ (•) is the measurement function based on the WSS yaw angle measurement described in (3.71),

• n a ∈ R ny+3 is a white noise n a ∼ N (0, R a t ).

The augmented measurement noise covariance matrix R a t is obtained from the block matrix concatenation of each of the different measurement noises:

R a t =      R t 0 0 0 R v t 0 0 0 σ 2 ψ      (3.75)
where R t was defined in (5.4), σ 2 ψ is the WSS yaw rate noise variance as defined in chapter 2, and R v t is given by

R v t =   σ 2 V 0 0 σ 2 cons   (3.76)
where σ 2 V is the WSS velocity noise variance as defined in chapter 2. In [START_REF] Shin | Estimation techniques for low-cost inertial navigation[END_REF] the noise variance of the constraint-measurement v y = 0 is related with some possible small misalignment angle ϑ between the IMU and the vehicle frame. The constraint-based noise variance is thus given by σ 2 cons ≈ v x t sin(ϑ), so for standard land vehicle applications it can be set to σ 2 cons ≈ 1m/s [God06].

Discrete state model

A discretization must be done to the GPS/WSS and GPS/INS/WSS continuous-time state models so they can be applied to the discrete EKF and UKF algorithms presented in table 3.1 and 3.3. The discrete state model can be written as:

X t = F t X t-1 + B t u t , (3.77)
where B t is the coefficient matrix used to shape white noise input so Q t = B t B T t and u t ∼ N (0, 1). Consider the state model is time invariant for the time interval T s over which the prediction is to be performed. As a consequence, the continuous-time state vectors (3.47) and (3.58) can be discretized by replacing derivatives by finite differences, e.g. ẋ(t) by x[(n+1)Ts]-x[(nTs)] Ts . The closed form simple solution for the discrete transition matrix is given by

F t = e Ts Ft .
(3.78) 92CHAPTER

HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

Furthermore, assuming that T s Ft ∞ is small enough the exponential can be approximated by e Ts Ft = I nx + T s Ft , so that the final discrete transition matrix is expressed as [START_REF] Nordlund | Sequential Monte Carlo Filters and Integrated Navigation[END_REF]:

F t = I nx + T s Ft (3.79)
Similarly, the covariance of the discrete noise sequence can be written as [START_REF] Nordlund | Sequential Monte Carlo Filters and Integrated Navigation[END_REF]:

E B t u t u T t B T t = T s I nx + T s 2 Ft Bt Qt Bt I nx + T s 2 F T t + T 3 s 12 Ft Bt Qt BT t F T t , (3.80)
where for practical reasons the last term is generally neglected, yielding the following explicit solution for B t ,

B t = t+Ts t I nx + (t + T s -τ ) Ft Bt dτ = T s I nx + T s 2 Ft Bt . (3.81)

Results

The two proposed hybrid approaches, GPS/WSS and GPS/INS/WSS, are tested using the EKF and UKF. Parameters tunning the UKF are listed in table 3.4 according to values used in [START_REF] Shin | Estimation techniques for low-cost inertial navigation[END_REF].

UKF parameters

k = 0 α = 0.001 β = 2
Table 3.4: UKF scaling parameters.

Performances are compared in terms of 2D horizontal error (HE). We recall here that the horizontal position error can be defined as functions of the latitude δλ and longitude δφ errors by

HE = (Rδλ) 2 + (Rδϕ cos (λ)) 2 , (3.82)
where R stands for the Earth's radius. The proposed navigation strategies are first validated using simulated data and later tested on real data.

RESULTS

Synthetic data

Several simulations have been conducted to compare the different coupling techniques. All results presented in this section have been averaged over 25 Monte-Carlo runs (note that the vehicle trajectories differ from one Monte Carlo run to another). An example of the simulated vehicle dynamics corresponding to an acceleration variance of 2m/s 2 is depicted in Fig. 3.5(a). The pseudo-range measurement accuracy for the GPS signal is 12m (i.e. the standard deviation of n ρ i in (3.50) equals 12m). The WSS parameters have been adjusted according to standard vehicles, i.e. L = 1.8m and R rr = R rl = 0.23m. INS errors have been simulated according to a Honeywell HG1700 tactical IMU [START_REF] Gao | Development of precise GPS/INS/Wheel Speed Sensor/yaw rate sensor integrated vehicular positioning system[END_REF] to outline the IMU contribution to the hybrid system. Its performances will be contrasted in the next section to a real low cost IMU. The spectral noise densities for the HG1700 [Pet03] inertial sensors are summarized in table 3.5. The standard deviations for the WSS noise components (σ V , σ ψ , σ rr , σ rl ) are summarized in the table 3.6 according to standard values given in [START_REF] Carlson | Practical position and yaw rate estimation with GPS and differential wheelspeeds[END_REF].

Gyroscope Accelerometers

σ g σ bg σ x a σ x ba σ y a σ y ba [ • /s/ √ Hz] [ • /s/ √ Hz] [m/s 2 / √ Hz] [m/s 2 / √ Hz] [m/s 2 / √ Hz] [m/s 2 / √ Hz]
1, 5.10 -3 1, 8.10 -4 4, 1.10 -3 4.10 -6 4, 1.10 -3 1, 1.10 -5

Table 3.5: HG1700 2D sensor noise parameters.

WSS

σ V σ ψ σ rr = σ rl [m/s] [ • ] [m] 0.04 3 10 -5
Table 3.6: WSS sensor noise parameters.

A scenario with full GPS visibility is first presented. The GPS receiver is assumed to view 7 satellites in line of sight (LOS). Figure 3.5(a) shows an example of real trajectory followed by the vehicle and the estimated position obtained when using WSS only (WSS trajectory). The estimated trajectory resulting from the coupling approach GPS/WSS (Filtered trajectory) is also depicted. The coupling between GPS and WSS is clearly necessary for this example. For further comparisons, a more realistic urban scenario is presented with partial GPS availability.

The idea is to observe the behavior of the different coupling techniques with a reduced number of visible GPS satellites. For this, we have simulated a 30s time interval during which only two GPS satellites are visible, followed by a 30s GPS outage (both indicated by vertical black lines in Fig. 3.7). 3.7 shows that the asymptotic horizontal error PCRBs corresponding to the GPS/WSS coupling tend to diverge over time when working with a small number of satellites, contrarily to the GPS/INS/WSS hybridization. Note that there is a significant gain when passing from 5 to 6 satellites in the GPS/WSS case. Indeed, the sixth satellite has a very good visibility (in terms of Dilution of Precision) for this GPS constellation. 3.7: Asymptotic Horizontal PCRBs.

RESULTS

95

RESULTS

Tests done on real data

The previous section was aimed at studying the theoretical performances and relevance of the different integration strategies. The PCRB was contrasted to the theoretical results to conclude on the correctness of the proposed system. An extension of these performances to real data is now analyzed.

In particular, tests will be done under challenging scenarios including urban canyon phenomena.

The experimental data was obtained from a test field campaign carried out in ISAE campus,Toulouse (France). This trajectory was specially chosen because it presented open sky areas and a semi-urban scenario phenomena. By semi-urban scenarios we understand that the measurement error sources (multipath, signal attenuation, partial or total GPS outages, etc...) are not as severe, in magnitude and duration, as for highly corrupted urban scenarios. The idea is to analyze the system behavior in conditions where the integrity of the GPS measurements is not assured (i.e. they do not behave according to their theoretical error model). The reference navigation solution was obtained from a synchronized position attitude navigation (SPAN) system composed by a Novatel receiver, with differential GPS (DGPS) approach, and a high-accuracy IMAR inertial unit. The reference trajectory is shown on white in Fig. 3.8. The followed path is 1.7 km long and it was done in 3.5 minutes. Odometric data was obtained from the already on board ABS WSSs. Data was recuperated through the on-board diagnostic (OBD) port using an ELM327 device as the one showed in chapter 2. The velocity frames were produced and processed at 50Hz. The inertial data was provided by a typical low cost MEMS-based IMU: Xsens MT9 [START_REF] Haykin | Motion Tracker Technical Documentation MTx-B[END_REF]. This unit works with analog components whose individual price is around 15USD. A 2 axis accelerometer and 1 axis gyro were used. Their frequency is 100 Hz. No significant loss of information was found when sub-sampling the INS data at 50Hz, so in this way the frequency of both dead reckoning systems was propagated at the same rate.

The correct modeling of the DR sensor errors presents a key point in the integration performance.

An underestimation of the sensor error parameters would prevent the GPS measurements from correctly compensating the sensor biases (because not enough importance would be given to pseudorange measurements). On the other hand, an overestimation of the error parameters would decrease the relevant contribution of the dead reckoning sensors to filter possible instantaneous anomalies in the GPS measurements. The next section describes the different methods employed in this thesis for the characterization of the DR sensors.

Sensor characterization

According to the error models given in chapter 2 for the WSSs and INS outputs, the following noise parameters need to be characterized:

• the wheel radius noises: σ 2 rl and σ 2 rr

• the noise of the WSS-derived yaw rate: σ 2 ψ

• the noise of the WSS-derived velocity:

σ 2 V          ∈ WSS measurements
• accelerometer noise: (σ x a ) 2 and (σ y a ) 2

• noise for the accelerometers' bias evolution:

(σ x ba ) 2 and (σ y ba ) 2    ∈ INS measurements
Indeed by WSS measurements we mean the velocity and yaw rate measurements derived from the WSS outputs.

Noise parameters for WSS measurements

Typical procedures to model the sensor errors include measuring the output of the static or idle sensor. In this way, no external time varying forces are considered to be significantly affecting the sensor output. The measured data is thus only representing the different error sources affecting the sensor. However, such characterization is not possible for the ABS-based WSSs. Wheel speed sensors (WSSs) used in this research always output zero speed when the vehicle is not moving. Hence no sensor characterization can be performed in usual static or idle conditions. A sub-optimal technique is thus proposed. An in-motion characterization is done where external factors possibly contributing to sensors' errors are minimized. The reference solution provided by the Novatel-IMAR-SPAN solution can provide velocity accuracy at a millimeter per second level and yaw-rate accuracy at millidegree per second. Considering that both WSS velocity and yaw-rate measurements usually present lower accuracies in at least one or two orders of magnitude, is suitable to consider the Novatel-IMAR-SPAN system as the navigation reference. To avoid sensing for errors in the nominal wheel radii or slippage effects, data collection for the WSS measurement characterization is done at slow velocities and for a straight trajectory. Results are presented in Fig. 3.10 for the velocity measurement and in Fig. are summarized in table 3.8. Due to the sub-optimal characterization technique, conservative noise parameters are to be considered herein.

On the other hand, it was out the scope of this thesis to test for real changes in the wheel radius under different conditions. According to [START_REF] Stephen | Development of a multi-sensor GNSS based vehicle navigation system[END_REF], a logical noise variance accounting for variations in the wheel radius can be set to σ rl = σ rr = 10 -5 m.

WSS σ

V σ ψ σ rr = σ rl [m/s] [ • ] [m] 0.1 3 10 -5
Table 3.8: WSS simulated sensor noise parameters.

RESULTS

Noise parameters for IMU sensors

To evaluate the noise terms affecting the INS sensors, data was collected at the output of the idle sensors during several hours. Under static conditions, the IMU measures the Earth's rotation rate and the force of gravity. However, as these two quantities remain constant during the data collection period, they will not interfere with the Allan variance computation (see chapter 2 for more details).

Moreover, for very low cost gyros as the one employed in this thesis, the Earth's rotation rate is usually masked by the sensor noise.

The Allan variance curves for the used Xsens MT9 accelerometers and gyroscope are respectively shown in Fig. 3.12 and Fig. 3.13. Evaluating the curves according to the principle explained in section 2.3 yields the sensor noise spectral densities presented in table 3.9. The used sensor bandwidth is 50Hz.

Comparing table 3.9 with table 3.5, it becomes evident that performances to be obtained with the real Xsens MT9 IMU will be poorer than the ones obtained in the previous section with a simulated HG1700-like IMU. In this way, the relevance of the contribution of the INS to the hybrid system will be analyzed in terms of the IMU quality.

[ • /s/ √ Hz] [ • /s/ √ Hz] [m/s 2 / √ Hz] [m/s 2 / √ Hz] [m/s 2 / √ Hz] [m/s 2 / √ Hz]
1.10 -1 8.10 -2 2.10 -3 3.10 -4 2.10 -3 4.10 -4

Table 3.9: Xsens MT9 2D sensor noise parameters.

Navigation Performance

A first study is done for a standard stand-alone GPS strategy using an EKF. Fig. 3.14 shows the horizontal 2D positioning error and Fig. 3.15 the number of visible satellites during the test field

RESULTS

campaign. A first phenomenon is observed during the first 50 seconds: a constant bias of approximately 10 meters is affecting the error solution. This undesired effect can be certainly expressed by the presence of erroneous pseudorange signals during the initialization period. Indeed, the starting point was surrounded by tall buildings, so it is highly possible that a constant multipath signal was tracked during the first seconds introducing a constant bias in the position solution. Two other interesting phenomena, highlighted by magenta circles, will be analyzed herein. Around time instant t = 100s and t = 135s important errors have a sudden impact on the navigation accuracy. From Fig. 3.16 it can be observed that these time intervals correspond to zones with tall buildings and/or with deep foliage. At the same time Fig. 3.15 shows that no significant reduction over the visible GPS constellation can justify the presence of important errors. Multipath interference, signal attenuation and cross-correlation errors will certainly be affecting the received GPS signals and thus the pseudorange measurements. As already explained in chapter 1, to distinguish the contribution of each individual phenomena to the signal's error presents major difficulties. Hence, all these effects degrading the signal quality will be summarized hereafter under the name of "measurement errors".

In the following, the contribution of real DR sensors to cope with all the above mentioned phenomena will be analyzed. in an ascending order and taking the median of the error set. Similarly, the 95% bound represents the lower threshold under which the 95% of the errors are found. It is important to understant that if just the means were compared for the different strategies, the presence or absence of high error peaks would get masked within the mean computation. In this way, for instance, the 95% bound provides useful information on the presence of error peaks. To facilitate the comparison, error values (in terms of mean and 95% bound) corresponding to the best performing strategy are highlited in green.

Similarly to results obtained from simulated data, no significant difference is observed between the UKF and EKF performances. Moreover, Table 3.10 shows that for the ensemble of the trajectory, the

RESULTS

stand alone GPS and GPS/WSS or GPS/INS/WSS approach present similar results. For scenarios assuming a good visibility over the GPS satellites, DR sensors do not really provide a meaningful contribution to the navigation solution. Indeed, in the presence of good quality GPS signals, the lowcost/low-accuracy DR navigation solution is not assigned an important "weight" during the filtering process. However, the DR contribution is put in evidence under degraded GPS reception conditions as summarized in Table 3.11. As a similar behavior was observed for the EKF and UKF, result are just given for the EKF. In the presence of "measurement errors" affecting the GPS signal, DR sensors will enable a smoothed or filtered position solution to be provided. Furthermore, results presented in Fig. 3.17 and the 95% bound in Table 3.11 show that the GPS/WSS strategy is better suited to smooth the sharp peak errors than the GPS/INS/WSS. It is recalled that in the GPS/WSS strategy the WSS fully determines the state vector while in the GPS/INS/WSS approach it is included in the measurement vector together with pseudoranges. It is also interesting to notice that from the comparison of the first 50 seconds in Fig. 3.14 and Fig. 3.17, it can be seen that DR sensors allow for a gradual compensation of the multipath interference during the initialization interval. 3.11: Horizontal error statistics in meters for the GPS degraded intervals using the EKF

• total GPS outage for 20s between t = 30s and t = 50s

• partial GPS outage for 15s between t = 120s and t = 135s

• total GPS outage for 20s between t = 150s and t = 170s were chosen for the partial GPS outage. Fig. 3.19 shows the different filtering and hybrid system performances for this scenario and Table 3.12 summarizes the results for the outage intervals. Once more the EKF and UKF applied to either the GPS/WSS or GPS/INS/WSS approach, attained

Conclusions

The principles of the different integration strategies were presented in this chapter. Two different filtering techniques were proposed for the navigation solution: the commonly-used Extended Kalman Filter (EKF) and the Unscented Kalman Filter (UKF). The UKF is aimed at overcoming the EKF limitations due to the linearization applied to the system model. In parallel, two dead-reckoning(DR) augmentations were proposed to the stand-alone GPS approach. Wheel speed sensors (based on differential odometry strategy) and a 2D IMU were studied as complementary navigation systems to the GPS. Two hybrid systems were proposed: GPS/WSS and GPS/INS/WSS.

Tests were done using both simulated and real data. Non linearities in the measurement model were tested during degraded GPS reception conditions. Non linearities in the state model were tested during GPS outages where DR error dynamics tend to be specially high. However, for none of these scenarios the UKF presented significant performance improvements over the EKF. Land vehicles are subjected to slow varying dynamics yielding similar performance for both techniques.

Tests done on simulated data validated the correct modeling of the hybrid strategies when contrasted to the Posterior Cramer Rao Bound (PCRB). The asymptotic PCRBs were also computed under varying numbers of LOS GPS satellites. The theoretical stability in terms of accuracy for the GPS/INS/WSS system appeared superior to GPS/WSS. Having more DR measurement sources enables the GPS/INS/WSS approach to be better adapted to work in urban areas affected by limited GPS availability.

The hybrid GPS/WSS and GPS/INS/WSS approaches showed important improvements over the stand alone GPS when tested for real data under degraded scenarios (i.e. in the presence of measurement errors affecting the GPS signal). Indeed, DR measurements helped to filter for abrupt changes (i.e. errors) in the position estimation induced by the presence of outliers. During GPS outages DR measurements ensured the continuity over the position estimation. While the WSSs presented highly interesting performances in real stand-alone mode, the relevance of the contribution of the INS to the WSSs navigation solution was shown to be subjected to the inertial sensors' quality. For simulated data, a good quality low cost IMU was considered. In this scenario, the GPS/INS/WSS strategy presented better performances. However, for the tested Xsens MT9 IMU, the coupling approach using just the WSSs (GPS/WSS) was preferred to the GPS/INS/WSS. Therefore, the GPS/INS/WSS strategy is discarded in the next chapters as a possible hybrid navigation solution. It is important to note, that this decision derives exclusively from the fact that the available low-cost IMU presented very poor performances. No generalization to the GPS/INS/WSS strategy has to be done.

In the previous chapter performances in terms of horizontal accuracy were given for the different hybrid strategies. However, no details on the reliability or consistency of the results were provided.

Indeed, no knowledge was given on how accurate the error models were adapted to the real error values. This is the purpose of the last two chapters of this manuscript. Conclusions on the performance of a navigation technique will not only be evaluated according to the final horizontal position error but also by considering the system capacity to provide a consistent confidence interval for the estimated position. In other words, a position solution must be given with a respective bound accounting for the accuracy of the estimation. In this chapter the robustness of Doppler measurements to urban canyon phenomena is going to be investigated. Indeed, this study is motivated by the interesting results shown in [START_REF] Kubrak | Vehicular navigation using a tight integration of Aided-GPS and low-cost MEMS sensors[END_REF] when including this type of measurement in the filtering model. Doppler measurements are going to be presented as a complement/ alternative to "noisy" pseudoranges. If the received pseudorange measurement does not appear to be consistent with its error model, Doppler measurements are going to be used for the navigation solution. In particular, Doppler measurements will be introduced in the system according to a "smooth correction criterion". Based on results presented in the previous chapter, the overall system performance will also be tested in the presence of DR navigation systems.

Different types of GPS measurements

Chapter 1 presented the models for different types of measurements that can be obtained after processing the received GPS signal. These measurements are:

• Pseudorange measurements,

• Phase measurements,

• Doppler measurements.

The properties related to each of these measurements are determined by the tracking loop used for their computation [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF]. A code synchronization is achieved by the delay lock loop (DLL) to obtain the pseudorange measurement. This type of tracking loop is very robust to dynamic stress but doesn't provide very accurate measurements (consider that synchronization is done over a code whose chip length is ≈ 300m). As the code tracking loop directly measures the signal time of arrival, the presence multipath will highly affect the tracking performance.

Phase measurements are obtained from a phase lock loop (PLL). PLLs replicate the exact phase and frequency of the incoming carrier signal to perform the carrier wipeoff function (i.e., convert the received modulated signal to its base bandwidth). Considering that the received signal is modulated in L1=1575.42MHz frequency, the phase synchronization is done over a signal wavelength of 0.19m providing very accurate cm-level measurements. Though this type of loop produces the most accurate velocity and positioning measurements, it is very sensitive to dynamic stress. Therefore it is not suitable for noisy scenarios such as urban environments or to perform the system initialization process. Moreover, phase measurement can be affected by cycle slips especially at low C/N 0 or high dynamics. Therefore the phase measurement cannot be exploited in degraded environments. Phase measurements are generally used for high accuracy applications in open-sky scenarios where steady state low or moderate dynamics can be guaranteed. It is important to note that even under very favorable conditions, pseudorange measurements are needed to initialize the phase measurement ambiguity term. Indeed, carrier smoothing techniques, as the Hatch filter [START_REF] Hatch | The synergism of GPS code and carrier measurements[END_REF] explained in section 4.2, are commonly used to propagate an initial pseudorange measurement by the phase measurement.

In this way, the initial ambiguity is resolved, and high precision positioning can be assessed.

The Doppler measurement can be still obtained in degraded scenarios thanks to the tolerance of the frequency lock loop (FLL) to high dynamic stress (a maximum FLL dynamic stress tolerance of 50Hz is assumed). FLLs perform the carrier wipeoff process by replicating the approximative carrier frequency, and they typically permit the phase to rotate with respect to the incoming carrier signal's phase. They provide less accurate velocity measurements than the PLL and absolute positioning cannot be achieved. Similarly to dead-reckoning systems explained in chapter 2, "relative" information given by Doppler measurements needs to be integrated to obtain the position solution. Hence, standalone Doppler navigation is not used because of its growing errors over time. However, it is an interesting measurement to be exploited in urban navigation because of the FLL robustness to noisy scenarios and to multipath interference [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF]. For instance, typical FLLs can track signals down to 23dBHz, for a noise bandwidth of 15Hz [START_REF] Van Dierendonck | GPS receivers[END_REF], while the PLL threshold is typically set around 30dBHz. The multipath "immunity" assumption is based on the fact that as the FLL performs the tracking over a high frequency signal (with a wavelength of approx. 20cm), the impact of multipath errors is highly bounded (contrary to the DLL case where the tracking is done over a code with chip length equal to 300m).

Considering the aforementioned characteristics for the three types of GPS measurements, a Dopplerbased smoothing approach will be discussed in this chapter. The idea is that degraded pseudorange measurements, mostly affected by urban canyon phenomena, may be complemented by more robust Doppler measurements. In particular, under good-visibility/low-interference GPS conditions, pseudorange measurements are considered to be reliable and they can be used for the position solution.

Otherwise, if one of these two conditions is not satisfied, Doppler measurements are exploited. One of the ideas explored in this chapter deals with the merging of pseudorange and Doppler measurements.

The principle for this merging is somehow linked to typically used carrier smoothing techniques briefly presented in the next section.

Carrier smoothing

A common processing technique that combines the absolute pseudorange measurement with the highly accurate carrier phase measurement is referred to carrier smoothing. A well known carrier smoothing algorithm is based on the Hatch filter [START_REF] Hatch | The synergism of GPS code and carrier measurements[END_REF] where the filter measurements are given by (neglecting the clock bias and noise term):

ρt,i = β t,i ρ t,i + (1 -β t,i ) (ρ t-1,i + (φ t,i -φ t-1,i )) , (4.1)
where

• ρt,i is the smoothed pseudorange at time instant t for the ith satellite,

• β t,i is the ith satellite weight term which varies from 1 to ≥ 0 by small increments (≤ 0.1),

• φ t,i is the phase measurement at time instant t for the ith satellite (1.10),

• (φ t,i -φ t-1,i
) is the range difference between two consecutive carrier phase measurements,

• ρ t,i is the ith measured pseudorange at time instant t (3.50).

The principle of this filter is the following: the pseudoranges are merged with the carrier phase measurements in a recursive filter that progressively increases the weight of the carrier phase while decreasing the weight of the pseudo-range measurement. The goal of this method is to smooth the pure pseudorange by the carrier phase measurement and eventually provide a smoothed range measurement that it is mostly obtained from the accurate carrier phase measurement. The time of transition between a highly dominated pseudorange and a highly dominated phase measurement is given by the filter smoothing time. This convergence time is set to be superior to multipath correlation time. Apart from the fact that in urban scenarios the user cannot rely on the continuity of the phase measurement, another limitation of this algorithm is the danger that the smoothed pseudorange ρ is initialized by a biased pseudorange ρ (for ex. with multipath interference). In this case, an erroneous measurement will be delivered. This error can be somewhat limited by keeping the weight β t at a value greater than zero so once the pseudorange is no longer affected by a multipath the initial introduced error can be eventually compensated.

Doppler smoothing

Based on the carrier smoothing technique we propose to use a new type of measurement for the navigation solution especially suited to urban scenarios. Contrary to the expression given in (4.1), we propose to use Doppler measurement to smooth pseudoranges. The Doppler measurement needs to be expressed as a pseudorange rate to enable one to merge pseudorange and Doppler measurements.

We recall that the Doppler effect is caused by the relative motion of the transmitting satellite with respect to the receiver. The expression given in chapter 1 for the Doppler measurement D t,i (1.8) is:

D t,i = - L1 c v e t -v S t,i T • l LOS i,t + d t + n D t,i . (4.2)
where

• L1 = 1575.42M Hz and c is the speed of light,

• v e is the receiver velocity vector in geodetic ECEF coordinates,

• v S i is the satellite velocity vector in geodetic ECEF coordinates obtained from ephemeris,

• l LOS i is the unit line of sight vector between the i th satellite and the receiver,

• d t is the receiver clock drift, • n D i,t ∼ N (0, σ 2 D t,i
) and σ 2 D t,i is mainly determined by (1.12).

This equation can be easily modified to achieve a range rate equation. Considering that the range measurement is given by:

r t,i = x t -x s t,i 2
+ y t -y s t,i 2
+ z t -z s t,i 2 , (4.3)
the range rate ṙ can be directly obtained by calculating the partial derivatives of (4.3) with respect to the vehicle position (x, y, z)

ṙt,i = ( ẋt -ẋs t,i )(x t -x s t,i ) + ( ẏt -ẏs t,i )(y t -y s t,i ) + ( żt -żs t,i )(z t -z s t,i ) x t -x s t,i 2
+ y t -y s t,i 2
+ z t -z s t,i 2 , (4.4)
where ȧ denotes the changing rate (i.e. the velocity) of a. Hence, comparing (4.2) with (4.4) we have that:

ṙt,i = v e t -v S i,t T • l LOS t,i , (4.5)
so the pseudorange rate derived from the Doppler measurement is written as:

ρ = - c L1 D t,i = v e t -v S t,i T • l LOS t,i + d t - c L1 n D t,i . (4.6)
The filter measurement obtained from the Doppler smoothing can now be expressed as

ρt,i = α t,i ρ t,i + (1 -α t,i ) ρt-1,i - c L1 D t,i ∆t (4.7)
where

• ρt,i is the smoothed pseudorange at time instant t for the ith satellite,

• α t,i is the ith satellite weight

• ∆t is the sampling time,

We will denote n ρt,i the noise corresponding to the ith smoothed pseudorange given in (4.7), so that n ρt,i ∼ N (0, σ 2 ρt,i ). The measurement noise variance σ 2 ρt,i is obtained as follows:

σ 2 ρt,i = α 2 t,i σ 2 t,i + (1 -α t,i ) 2 σ 2 ρt-1,i + c L1 σ D t,i ∆t 2 (4.8)
where σ 2 t,i is the noise variance associated to the ith pseudorange measurement at time instant t. Contrary to section 4.2, eq. (4.7) proposes an approach where the weights will not be time-correlated. Indeed, their value will be constantly estimated according to the pseudorange measurement reliability.

Note that for the initialization of ρt,i α t,i is forced to 1.

Measurement reliability

Reliability is described in [START_REF] Petovello | Real-Time Integration of a Tactical-Grade IMU and GPS for High-Accuracy Positioning and Navigation[END_REF] as the ability to identify and reject non-modeled errors in the measurements. These erroneous measurements, called outliers, will bias the navigation solution and thus it is important to detect any anomaly in the observation and provide an adequate processing strategy. The actual identification of outliers is performed via testing of the measurement residuals in

RAIM+FDE strategy

Integrity is generally defined as a measure of the trust that can be placed in the correctness of the information supplied by a system. It also includes the availability of the system to provide timely and valid warnings to the user when it is inadvisable to use the system for navigation. This study particularly focuses on the first of these two characteristics describing integrity. An error bound is a measure of trust on the accuracy of the estimation being done. Among the different existing integrity procedures, the snapshot RAIM is the best known and widely implemented approach [C + 06a]. The FDE is a complement to this system where not only the presence of an anomaly is detected but its source (i.e. the corresponding satellite) is identified and further excluded from the navigation solution. The weighted least square (LS) estimator of the state vector at time instant t, derived from the linearized measurement model is classically defined as

X t = X t-1 + ∆ X t , (4.9)
with

∆ X = (H T t R -1 t H t ) -1 R -1 t H T t ∆Y t , ∆Y t = Y t -h( X t-1 ),
where H t is the linearized measurement matrix around X t-1 , R t is the measurement noise covariance matrix and X t-1 is the initial value of the state vector. The vector of LS residuals is thus defined as:

w t = Y t -h( X t ) = ∆Y t -h( X t ) -h( X t-1 ) = I ny -H t (H T t R -1 t H t ) -1 R -1 t H Tt ∆Y t .
Different test statistics can be used with the snapshot RAIM strategy to detect an anomaly. The range comparison method [START_REF] Lee | Analysis of range and position comparison methods as a means to provide gps integrity in the user receiver[END_REF], the least square (LS) residual method [START_REF] Parkinson | Autonomous GPS integrity monitoring using the pseudorange residual[END_REF] and the parity method [START_REF] Sturza | Navigation system integrity monitoring using redundant measurements[END_REF] are three of the main snapshot RAIM algorithms that use different test statistics. Due to its low computational cost, the LS residual method is going to be used in this paper. The corresponding test statistics is defined as:

SSE = w t R -1 t w T t = w t 2 . (4.10)
where SSE stands for sum of the squared errors. The test statistics SSE follows a χ 2 distribution with n y -4 degrees of freedom [START_REF] Martini | An integrity monitoring technique for multiple failures detection[END_REF]. Therefore a minimum of 5 satellites is needed to perform the integrity test. However its detection capability presents many drawbacks. For instance, as the state vector used to compute the residuals was initially computed using the received "defective" measurements, an assimilation of the error might already be present in X t . Especially in the occurrence of large or multiple biases, a compensation might be done within the residuals so the abnormal situation is finally not detectable. This occurs because when computing the residual term (4.10), a projection of ∆Y t is done on the null space of H T t [START_REF] Martini | An integrity monitoring technique for multiple failures detection[END_REF]. If the combination of measurement errors is such that the angle between ∆Y t and this subspace is big enough, many information will be lost in the projection. The norm of the residual vector w will be small and the error detection will fail.

The exclusion procedure is performed after an error detection has been achieved. The fundamental assumption used by the snapshot RAIM+FDE approach is that just one outlier can be present. A total of n y subsets each containing n y -1 satellites is built. For this purpose a minimum of 6 satellites must be available. A fault detection condition is performed over the n y subsets as in the detection step. Just the subset excluding the outlier will be under a given fault detection threshold. Therefore, this subset will be kept to compute the final position estimation. The covariance matrix P t associated to the estimate state vector X t is given by:

P t = (H T t R -1 t H t ) -1 . (4.11)
Considering that a correct detection and exclusion was achieved for the received measurements, matrix P t should be accurately modeling errors in the estimated navigation solution. Though the snapshot RAIM+FDE approach is widely accepted for civil aviation applications, its principles of functioning and specifications are not appropriate for land vehicle navigation. Due to either the presence of multiple outliers or to a limited GPS satellite constellation, the snapshot RAIM+FDE system is often unavailable in urban scenarios.

In this chapter signal reliability in urban scenarios is addressed not by excluding erroneous measurements but by replacing them by other types of measurements that may be less affected by urban phenomena. For this approach a minimum number of satellites is not required because each measurement is individually tested for anomaly presence using EKF innovations. The idea is that depending on the consistency between the received measurement and its associated error model, different types of measurements are going to be exploited. In particular, the smoothed pseudorange based on Doppler measurements (4.7) is going to be used.

Smooth correction function

Reliability testing results will determine the type of measurement to be used for the navigation filter.

We recall the expression for the EKF innovation vector I t (table 3.1):

I t = Y t -Ŷt|t-1 , (4.12)
where

• I t = (I t,1 , . . . , I t,ny ),
• Y t is the measurement vector Y t = (Y t,1 , . . . , Y t,ny ),

• Ŷt|t-1 = h( Xt|t-1 ) is the vector containing the pseudorange measurements predicted from the propagated state vector.

Assuming the state and measurement noises follow a zero mean Gaussian distribution and that their covariance matrices can be accurately characterized, the innovation covariance S t matrix is given by:

S t = H t P t|t-1 H T t + R t , (4.13)
where

• H t is the first order Taylor series expansion of the non linear measurement function h t ,

• P t|t-1 is the a priori state covariance matrix,

• R t is the noise covariance matrix associated to Y t .

The innovations associated to different time instants are assumed to be independent. The variance of the ith innovation at time instant t (denoted as s 2 t,i ) is given by the ith element of the diagonal of S t . As a consequence, in nominal situations (i.e. when measurement models are consistent with the observations) the ith innovation error is distributed according to a zero-mean Gaussian distribution with variance s 2 t,i , i.e. with pdf

p(I t,i ) = 1 √ 2πs t,i exp - 1 2s 2 t,i I 2 t,i . (4.14)

MEASUREMENT RELIABILITY 123

The EKF minimizes the mean square error (MSE) of the state vector considering that the received measurements have a Gaussian distribution with known parameters. If such a hypothesis is not fulfilled, the convergence and stability of the filter are no longer guaranteed. Hence, it is very important to detect any inconsistency in the filter models.

An indicator of received pseudorange quality will be obtained by comparing the values obtained in (4.12) with their assumed statistical distribution given by (4.13). It must be pointed out that any lack of consistency in the innovation sequence will be attributed to the uncorrect modeling of the measurement noise distribution given by R t . Neither the state model nor the linearized measurement transition matrix will be considered to introduce significant errors. We recall that noise variances σ 2 t,i included in R t are obtained as a function of the C/N 0 parameter and EGNOS corrections as explained in chapter 1, section 1.2.2. In particular, in the absence of temporal correlation in the pseudorange measurements and among the receiver channels, R t is a diagonal matrix fully described by σ 2 t,i . In open-sky scenarios, it can be reasonably considered that pseudoranges follow a zero mean Gaussian distribution where the C/N 0 parameter provides reliable information on the signal quality.

However, due to the particular properties of urban canyon phenomena, "measurement errors" are not fully captured by the C/N 0 parameter. In this case, the noise terms are usually underestimated and the pseudorange measurement model determined by N (0, σ 2 t,i ) is no longer valid. The test statistic usually employed to detect the presence of an anomaly in the ith received measurement is given by:

T t,i = |I t,i | s t,i (4.15)
where |a| denotes the absolute value of a. The test is normally compared to an established threshold • T t,i < ζ 1 the measurement model described by σ 2 t,i is considered correct and the ith pseudorange measurement is used to update the state vector.

ζ. If T t,i <
• T t,i > ζ 2 the received pseudorange is considered unreliable with respect to its noise model and possibly affected by multipath, cross-correlation errors, etc. The Doppler measurement is used to update the state vector. For a typical one-threshold-based anomaly detection is done for 3 < ζ < 6 [START_REF] Bullock | Integration of GPS with other sensors and network assistance[END_REF]. We propose herein to set ζ 1 = 2 and ζ 2 = 4. This choice has been motivated by results obtained for several sets of real data as presented in [START_REF] Spangenberg | Urban navigation system for automative applications using HSGPS, inertial and wheel speed sensors[END_REF]. If the innovation is smaller than twice its associated variance, which accounts for the 95% of the cases if the noise model is correct, the pseudorange is exploited. Otherwise, the pseudorange is totally ignored if the innovation is greater than four times its associated variance, which corresponds to a fault detection probability of 0.06%. The smooth correction function used for the weight computation is given by:

α t,i = - 1 ζ 2 -ζ 1 (T t,i -ζ 1 ) + 1 ∀ ζ 1 < T t,i < ζ 2 . (4.16)
Indeed, it can be also extended to T t,i < ζ 1 where the weight is constantly set to α t,i = 1 (Fig. 4.2). Therefore, ρt,i = ρ t,i and the pseudorange measurement is used for the filter update. However, the extension of the function to T t,i > ζ 2 , implying that the weight is constantly set to α t,i = 0 for this interval, is not considered. If α t,i = 0 the used measurement, in the absence of bias and noise terms, would be given by ρt,i = ρt-1,i -c L1 D t,i ∆t . This would mean that the stand-alone Doppler measurement is used to estimate position and not velocity parameters. Therefore integrated errors would unnecessarily be propagated. Moreover, we will need to assume that a reliable previous ρt-1,i exists, which might not always be the case. In conclusion, if T t,i > ζ 2 Doppler measurements are used to estimate the vehicle velocity according to the model given in eq. (4.2). To ensure the accuracy and reliability of the chosen measurements some constraints must be considered:

1. Stand-alone Doppler measurements are used to estimate velocity. If tracked for long periods, integration errors will get accumulated and divergent position estimations will be obtained. We propose to allow the use of consecutive Doppler measurements for the same satellite for a period not exceeding the 5 seconds. Otherwise the system will switch to the pseudorange measurement providing absolute positioning.

2. For the smoothed pseudorange expression (4.7) to be valid for ζ 1 < T t,i < ζ 2 , a previous reliable pseudorange ρ t-1,i or smoothed pseudorange ρt-1,i is needed to propagate the Doppler measurement. If the measurement used for the ith satellite in t -1 was a Doppler measurement, eq. (4.7) cannot be used. In this case the system will switch to the Doppler measurement (if this situation is not in contradiction with the previously described constraint) .

P t = (I nx -K t H t ) P t|t-1 , (4.17)
where the expressions for H t , K t and P t|t-1 can be easily derived from the filter state-space models presented in the following sections.

Ft =      0 2×2 A 0 2×2 0 2×2 0 2×2 0 2×2 0 2×2 0 2×2 F GP S t      , A =   1 R λ +ht 0 0 1 (Rϕ+ht) cos(λ)   , Qt = diag(0 2×2 , σ 2 n , σ 2 e , σ 2 b , σ 2 d ). ( 4

.19)

where A was described in (2.12), and F GP S t and (σ 2 b , σ 2 d ) were defined in (3.46). The measurement model (3.2) is redefined at every time instant depending on the type of measurement being used for each of the n y visible satellite. The measurement vector is defined as:

Y t = (Y t,1 , . . . , Y t,ny ) T , (4.20)
where Y t,i , i = 1 . . . n y , can be either a pseudorange ρ t,i , smoothed pseudorange ρt,i or Doppler measurement D t,i . According to this, the measurement function h t will be constructed as the concatenation of the individual measurement functions describing the exploited measurements:

h t =      h t,1 . . . h t,ny      , where h t,i =    q t,i • g(λ t , ϕ t ) + b t if Y t,i = ρ t,i or Y t,i = ρt,i p t,i • k(v n,t , v e,t ) -L1 c d t if Y t,i = D t,i (4.21)
where q t (•) was given in (3.54), and g t (•) is obtained from (3.51) considering that the position estimates are given by the state vector variables (λ t , ϕ t ). Similarly to the pseudorange model, a frame transformation must be done so the Doppler measurement can be written as a function of the state vector parameters. The function for the Doppler measurement p t is given by:

p t,i (v e t ) = - L1 c v e t -v S t,i T • l LOS i,t (4.22)
where v e t = (v x,t , v y,t , v z,t ) T is the vehicle velocity in the ECEF frame. Hence, k(v n,t , v e,t ) denotes the frame transformation matrix from the NED to ECEF frame expressed as,

k(v n,t , v e,t ) =      v x,t v y,t v z,t      =      -sin(λ t ) cos(ϕ t ) -sin(ϕ t ) -sin(λ t ) sin(ϕ t ) cos(ϕ t ) cos(λ t ) 0        v n,t v e,t   . (4.23)

Enhanced hybrid filter using smooth correction criterion

The smooth correction criterion is going to be implemented using DR sensors to provide a final navigation system that is well-adapted to urban scenarios. DR data can provide a highly informative reference frame to achieve a reliability monitoring on the GNSS signals. The "autonomous integrity monitoring extrapolation" [START_REF] Diesel | GPS/IRS AIME: Certification for sole means and solution to RF interference[END_REF] and the "multiple solution separation" [START_REF] Brenner | Integrated GPS/inertial detection availability[END_REF] are examples of integrity schemes that use inertial information to improve fault detection and exclusion functions.

However they were developed according to civil aviation requirements and their functioning and assumptions are not directly suitable for land vehicle navigation.

Based on results obtained in chapter 3 the well performing wheel speed sensors (WSS) are going to be used. Details on the WSS hybrid system based on differential odometry were given in section where state vector parameters (δV t ,δψ t ) represent the along-track velocity and yaw angle error associated to the corresponding WSS measurements.

A further enhancement is proposed in this chapter to deal with inaccurate yaw angle measurements provided by differential odometry techniques. A 1 axis gyro specially designed for land vehicle navigation is included in the system. In this way, WSSs will only provide along-track velocity estimations (i.e., estimations of V t ) while the heading information ψ t will be totally determined by the gyro outputs. ),

• δV t is the WSS along track velocity error and (δR rl t , δR rr t ) are the wheel radius errors,

• δψ t is the gyro yaw angle error and b g,t is the gyro bias,

• ( ḃt , ḋt ) are the receiver clock bias and clock drift.

The state transition matrix and noise matrices are as follows:

F GPS/WSS/Yaw t =   F W SS/Y aw t 0 0 F GP S t   , (4.27)
F WSS/Yaw t =                  0 0 cos(ψ) R λ +h V -sin(ψ) R λ +h 0 0 0 V tan(λ) 2 Rϕ+h 0 sin(ψ) (Rϕ+h) cos(λ) V cos(ψ) Rϕ+h 0 0 0 0 0 0 0 ω rl 2 ω rr 2 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0                  , BGPS/

Details on bound computation

To study the consistency of the proposed strategies, the horizontal position error (3.82) will be presented together with its respective bound. The horizontal position bound will be obtained from the updated state covariance matrix P t of the EKF. Consider the matrix P H,t containing only the variance and covariance parameters for the horizontal position estimates,

P H,t =   σ 2 λ cov(λ, ϕ) cov(ϕ, λ) σ 2 ϕ   (4.29)
where cov(a, b) denotes the covariance between a and b. As the horizontal error is described in the north(n)-east(e) frame, a suitable frame transformation must be applied to P H,t as follows

PH,t =   σ 2 λ 1 (R λ +h) 2 cov(λ, ϕ) 1 (R λ +h)(Rϕ+h) cos(λ) cov(ϕ, λ) 1 (R λ +h)(Rϕ+h) cos(λ) σ 2 ϕ 1 (Rϕ+h) 2 cos(λ) 2   =   σ 2 n cov(n, e) cov(e, n) σ 2 e   (4.30)
where PH,t is the covariance error matrix in meters for the horizontal north(n)-east(e) plane and all the frame transformation parameters were already described in (2.31). Lets consider X H = (δn, δe) T the vector containing the n -e position errors (δn, δe) where the time index is not included for notational simplicity. The general expression of X H pdf is given by:

p(X H ) = 1 2π det( PH ) e -1 2 X T H P -1 H X H (4.31)
where det(A) denotes the determinant of matrix A. The vector X H is a bi-dimensional random variable, X H ∼ N 0, PH . As PH is a positive definite matrix, it can be diagonalized and its eigenvalues are all positive. In particular, we can find an orthonormal basis B = ( ν 1 , ν 2 ) composed of eigenvectors ( ν 1 , ν 2 ) corresponding to the eigenvalues (κ 2 1 , κ 2 2 ) of PH . Hence, PH can be expressed as

PH = GDG T (4.32)
where

• D = diag(κ 2 1 , κ 2
2 ) is the diagonal matrix whose elements are eigenvalues of PH ,

• G is the projection matrix whose columns are the eigenvectors ( ν 1 , ν 2 ). In particular, G is orthogonal so

G T = G -1 .
Let X ⊥ be the projection of X H in the orthonormal basis B, so X ⊥ = G T X H . X ⊥ is a 2-dimensional Gaussian vector with covariance matrix D so that X ⊥ ∼ N (0, D). It can be easily observed that as the off-diagonal terms of D equal zero, the components of X ⊥ are mutually independent. Recalling that the horizontal error is obtained as

HE = (δn) 2 + (δe) 2 = X ⊥ , (4.33)
it is of interest to introduce the normalized square magnitude of X ⊥ given by X ⊥

2 n = X T ⊥ D -1 X ⊥ . In the general case, X ⊥ 2 n follows a χ 2
2 pdf with 2 degrees of freedom if κ 2 1 and κ 2 2 are different from zero.

The density function presented in (4.31) defines a two-dimensional bell-shaped surface. Contours of constant density are obtained by setting the exponent (X T H P -1 H X H ) to a constant. Indeed, considering that (X

T H P -1 H X H ) = X ⊥ 2 n
, this contour will describe an ellipse whose semi-major axis is oriented along ν 1 and whose semi-minor axis is oriented along ν 2 . The analytical expression of the ellipse is:

X ⊥ (1) 2 a 2 κ 2 1 + X ⊥ (2) 2 a 2 κ 2 2 = 1 (4.34)
where X ⊥ = (X ⊥ (1), X ⊥ (2)) T and a ranges over positive values. The parameter a is defined as

p( X ⊥ 2 ≤ a 2 ) = 1 -P c. (4.35)
where P c represent the tolerated probability for the real error not to be contained within the elliptical contour. In the general case, a 2 is obtained by evaluating the inverse cumulative distribution of the χ 2 2 pdf for P c . However, accounting on usual characteristics presented by the navigation error pdf, a worst case bound is calculated herein. It was shown in [START_REF] Brenner | Navigation system with solution separation apparatus for detecting accuracy failures[END_REF] and [START_REF] Vanderwerf | Fde using multiple integrated GPS/Inertial Kalman filters in the presence of temporally and spatially correlated ionospheric errors[END_REF] that one of the eigenvalues

(κ 2 1 , κ 2
2 ) usually dominates over the other. This means that the contour ellipse defined in (4.34) will be strongly orientated in the direction of the dominant eigenvalue. Just as an example, and without any loss of generality, lets consider that the dominant eigenvalue is κ 2 1 , so κ 2 1 >> κ 2 2 . In this case, the expression in (4.34) becomes

X ⊥ (1) 2 a 2 κ 2 1 = 1 ⇒ X ⊥ (1) aκ 1 = 1 (4.36)
Recalling that X ⊥ (1) κ 1

∼ N (0, 1), the value of a is obtained from (4.35) and (4.36) as:

a = Q -1 1 - P c 2 , (4.37)
where Q -1 (•) represents the inverse cumulative distribution function of the Gaussian pdf N (0, 1). In [START_REF] Brenner | Navigation system with solution separation apparatus for detecting accuracy failures[END_REF] and [START_REF] Vanderwerf | Fde using multiple integrated GPS/Inertial Kalman filters in the presence of temporally and spatially correlated ionospheric errors[END_REF] authors propose to model the contour of constant density associated to (4.31) as an overbounding circle of radius aκ 1 (always considering κ 1 the dominant eigenvalue). Fig. 4.4

illustrates the principle. In a general way, the final 2D horizontal error bound B is given by the circle

X ⊥ (1) 2 r 2 + X ⊥ (2) 2 r 2 = 1, (4.38)
where the circle radius r is define as

r = max κ 2 1 , κ 2 2 a. (4.39)
and a is obtained from (4.37). This bound is going to be used in the next section to analyze the consistency of the positioning solution.

Results

Gyro noise parameters

The Allan variance curve for the MLX90609 Melexis gyro is shown in Fig. 4.7. According to the principle explained in section 2.3, the resultant sensor noise parameters are given in table 4.1. Considering the low accuracy of yaw angle measurement obtained from the WSS differential odometry approach (table 3.8), it is immediate to conclude on the significant contribution the additional gyro will have to the GPS/WSS/Yaw navigation solution.

Gyroscope

σ g σ bg [ • /s/ √ Hz] [ • /s/ √ Hz]
1.10 -3 9.10 -4

Pseudorange vs. Doppler analysis

In this section, the impact of urban phenomena on pseudorange and Doppler measurement errors is going to be studied. In particular, special attention is paid to the consistency between the errors and the measurement noise models. The analysis is going to be done for two different satellites that were visible during the test field campaign. pdf given in red.

From the comparison between the pseudorange and Doppler pdfs, it can be observed that Doppler measurements are generally better bounded by their associated error model. Fig. 4.9(c) presents a tailed histogram while the histogram in Fig. 4.9(d) is better bounded when compared to the reference Gaussian pdf. Indeed, for this high elevation satellite the Doppler error model seems to slightly overestimate the actual errors. Similar observations can be done for satellite 5. Indeed, as satellite 5 has a low elevation angle, it is usually more susceptible to multipath and cross-correlation errors than satellite 31. Fig. 4.10(c) shows that pseudorange noise parameters highly underestimate the real error. We recall that noise variances are computed from the signal's C/N 0 and EGNOS corrections (section 1.2.2). Therefore they don't take into account the presence of "measurement errors". If errors are better bounded for one measurement than for the other one, conclusions can be drawn on the robustness of the former measurement to urban phenomena. In this way, the overall results validate the correct bounding of Doppler noise models as well as the assumption on Doppler robustness to urban phenomena.

Navigation performance

Performances using the smooth correction criterion are first going to be presented for a stand-alone GPS filter and later for the hybrid GPS/DR approach. Based on results given in the previous chapter, an EKF is used. For the GPS filter the noise variances modeling the state parameters vn,t and ve,t are set to σ 2 n = σ 2 e = 2m/s. The bound computation described in section 4.7 is done for P c = 10 -5 . Results obtained with the standard EKF, the snapshot RAIM+FDE strategy and the proposed smooth correction criterion, using only GPS signals for the navigation, are depicted in Figs. 4.11, 4.12(a) and 4.13. The number of measurements replacing nominal pseudoranges for the smooth correction criterion are depicted in Fig. 4.14. errors highlighted in green correspond to time instants where the reliability approaches were active. For the snapshot RAIM+FDE strategy this means that the outlier measurement (i.e. pseudorange) was successfully excluded from the solution. For the smooth correction criterion it means that at least one satellite pseudorange measurement was considered not reliable and replaced by an alternative measurement (i.e. smoothed pseudorange or Doppler). According to all these figures the following observations can be done:

• The standard EKF implemented without the proposed smooth correction criterion presents significant errors and inconsistent solutions (that are not within the computed bounds).

• The widely applied least square snapshot RAIM+FDE solution presents several limitations.

The system is unavailable during 10% of the time inducing very important errors. For instance, three main errors may be outlined around the time instants t = 600s, t = 950s and t = 1200s.

In the first case, the RAIM+FDE does not performed satisfactorily because only 4 satellites are in LOS. For the other two errors, the exclusion function could not be performed because multiple outliers were present among the received GPS signals. In this case, though a problem was detected, none of the n y -1 satellite subsets could be considered bias-free, which prevented the use of the FDE function. Moreover, a bound that largely exceeds the real errors entails an excessively big confidence interval.

• The navigation filter using the smooth correction criterion in Fig. 4.13 presents highly interesting performances. If compared with Fig. 4.11, it can be seen that most of the inconsistent

Bottom figure: Doppler measurements

A further correction strategy is implemented to put in evidence the contribution of the innovative smoothed pseudorange. Consider that the correction criterion is no longer described by a two threshold approach, but by a sole cut-off threshold. In this way, either the pseudorange is reliable, either it must be replaced by its corresponding Doppler measurement. No gradual transition, represented by correctly described (i.e. contained) by the computed bound. Error values highlighted in green represent the best performances.

GPS filter

The smooth correction criterion is now applied to a hybrid GPS/DR-sensors approach to provide a complete navigation system especially adapted to urban scenarios. Based on the good performances observed in the previous chapter, wheel speed sensors (WSSs) are used. An additional 1 axis gyro specially adapted to vehicle navigation is tested to enhance the DR navigation performance. Fig.

Conclusions

This chapter presented an innovative technique to work on the reliability and integrity of GPS measurements in urban scenarios. The proposed smooth correction criterion was aimed at providing an intelligent strategy to exploit the positive aspects of the pseudorange and Doppler measurements. In particular, "relative " Doppler measurements were proposed as an alternative/complement to "absolute" pseudorange measurements. Performances were contrasted with the commonly used RAIM+FDE integrity strategy. As the reliability analysis is individually performed on each measurement for the smooth correction criterion , no satellite visibility constraints need to be met. In this way, the proposed integrity strategy did not only overcome the RAIM availability problem but it also enabled a low error, accurate bounded solution. The GPS/WSS/Yaw hybrid system using the smooth correction criterion appeared as an excellent choice when accurate and reliable position estimations need to be guaranteed under challenging scenarios.

rely on the robustness of alternative GPS measurements (hypothesis that presented some flaws in the previous chapter). Herein, only pseudorange measurements will be exploited. Indeed, in the presence of outliers an innovative algorithm will be deployed to compensate the error presence.

The navigation system is complemented by a two-step detection procedure that aims at classifying outliers according to their associated source of error. Two different situations will be considered in the presence of multipath. These situations correspond to the presence or absence of line of sight signal for the different GPS satellites. Therefore, two kinds of errors are potentially "corrupting" the pseudoranges, modeled as variance changes or mean value jumps in noise measurements. An original multiple model approach is proposed to detect, identify and correct these errors and provide a final consistent solution. The GPS/WSS/Yaw hybrid system which showed interesting navigation performance in the previous chapter will be also tested with the proposed multipath mitigation strategy. Indeed, we are going to present an approach where we will benefit from dead-reckoning sensor multipath immunity to enhance the outlier detection.

Overview on existing multipath mitigation techniques

Several methods can be found in the literature concerning multipath mitigation. Different configurations of antenna arrays are among the hardware solutions [START_REF] Counselman | Multipath-rejecting GPS antennas[END_REF], [START_REF] Williams | Multipath mitigation performance of planar GPS adaptive antenna arrays for precision landing ground stations[END_REF]. Working on the receiver correlator output to mitigate the impact of multipath during tracking is another well known approach [DFF92], [START_REF] Macabiau | Use of multicorrelator techniques for interference detection[END_REF]. However, all these strategies are characterized by their high complexity. In order to avoid these difficulties (and be hardware-independent), multipath mitigation can be performed on the data processing stage. Each "defective" pseudorange is considered to be affected by an error that represents the total contribution of all the multipath signals to the measurement computation.

Multipath detection and correction has already received some attention in the literature. Giremus Gaussian mixtures and the tracking was performed using particle filtering. However, the existing algorithms described above require to define a priori distributions for the NLOS error. This a priori knowledge is not easy to obtain in real urban scenarios. Moreover, the high computational cost of particle filters is a problem for land vehicle applications.

Multipath interference

As already explained in section 1.2.3, in the presence of multipath pseudoranges are affected by errors.

The most important errors occur when the multipath is in phase with the direct signal (i.e. φ1 = 0 • ) or in phase opposition (i.e. φ1 = 180 • ). It is important to note that the error magnitude depends on the correlation function shape and on the type of discriminator used. The error envelope presented in Fig. 5.1 shows the maximum errors attained in the GPS pseudorange calculation as a function of the multipath delay interfering the direct signal for a given receiver configuration. Any error between these two bounds, depicted by the shaded area, corresponds to -180 • < φ1 < 180 • . A coherent dot product delay lock loop (DLL) was used for the code tracking with 1/2-chip spacing between each correlator. The multipath was assumed to be attenuated by 3 dB after reflections (i.e. Ã1 = 0.5).

For a deeper insight on the GPS signal tracking process the reader is invited to consult [START_REF] Ward | Satellite signal acquisition, tracking, and data demodulation[END_REF]. In order to detect the presence of an anomaly in the GPS signal, the actual error distribution will be studied and contrasted to the nominal error distribution due to the sole presence of noise.

Therefore, several pseudorange samples will be needed. Though instantaneously a pseudorange might seem affected by a bias, within an observation window, this error can remain constant or vary over time. This phenomenon will depend on the characteristics of the received reflected signal at each time instant, that themselves depend on the navigation scenario and vehicle dynamics. This chapter proposes to study two kinds of errors: either the error remains approximately constant over the observation window and is modeled as a mean value jump [GTC07; HV06], either the error varies over time and it is modeled as a noise variance jump (in [MHG + 05] this phenomenon is studied for aerial navigation). To validate these two hypotheses Fig. Hereafter, the presence of multipath introducing a mean value jump will be referred to as NLOS interference, while the noise variance jump will correspond to the LOS interference. In NLOS situations, only a reflected signal due to multipath is received and tracked (a constant bias is present in the pseudorange measurement). In the LOS situation, the measurement is composed by the direct signal plus delayed reflections. Considering that the relative parameters (i.e. Ã, φ, τ ) vary over time, the value and sign of the introduced errors will be constantly changing and this situation will be modeled by a noise variance jump. Indeed, the actual mean value jump does not only represent a NLOS situation but also an LOS situation where the vehicle is not moving (i.e. the multipath amplitude, delay and phase are constant, so its final contribution to the direct received signal is also a mean value jump). However, as this error is finally considered under the hypothesis H 1 , the approach is still valid without loss of generality. The following section makes explicit the relation between these two types of errors and the received measurement model.

State-space models

State model

The state model described in section 4.5 is used herein. We recall that the state vector is given by: Ẋt = λt , φt , vn,t , ve,t , ḃt , ḋt

T GPS ∈ R 6 , (5.1)
For further details on the state transition and noise matrix please refer to section 4.5.

Measurement model

The pseudorange measurements associated to the ith satellite can be written as • b t is the receiver clock bias,

Y t,i = p s t,i -p t + b t + m t,i + (σ 2 t,i + r 2 t,i )w t,i , (5.2
• m t,i represents a possible mean value jump for the ith satellite,

• σ 2 t,i is the measurement noise variance in nominal conditions for the ith satellite,

• r 2 t,i represents a possible noise variance jump for the ith satellite,

• w t,i is a zero mean Gaussian variable such that w t,i ∼ N (0, 1).

The n y measurements are usually concatenated according to

Y t = h(X t ) + M t + R 1/2 t w t (5.3)
where

• Y t = (Y t,1 , ..., Y t,ny
) T is the pseudorange measurement vector,

• h(X t ) = p s t -p t + b t represents the measurement function including a frame transformation from the geodetic state coordinates to the rectangular coordinates as described in (3.51),

• M t = (m t,1 , ..., m t,ny ) T where m t,i = 0 if the ith satellite is affected by a mean jump and m t,i = 0 otherwise.

Assuming the pseudoranges are independent, the measurement noise covariance matrix R t is expressed as

R t =         σ 2 t,1 + r 2 t,1 0 • • • 0 0 . . . . . . . . . . . . 0 0 • • • 0 σ 2 t,ny + r 2 t,ny         . ( 5.4)
where r t,i = 0 if the ith satellite is affected by a noise variance jump and r t,i = 0 otherwise.

Multipath mitigation approach

Multi-hypothesis approach

Usual multipath mitigation schemes based on pseudorange measurements consider a binary system where the received signals are either bias-free or they are subjected to a multipath interference. In this chapter we propose to further develop the interference processing by introducing two different models for LOS and NLOS interferences. The assumption considering an error introduced by the presence of multipath is decomposed in order to identify the specific source of this error. Although a change in the noise variance (i.e. NLOS case) does not have such a strong impact on the positioning accuracy when compared to a mean value jump, it reveals to be a crucial factor when a precise bounding must be given for the final position solution. As a consequence, we propose a three hypothesis model to detect, identify and correct measurement errors due to multipath:

• H 0 : absence of error (only the direct signal is tracked and nominal σ 2 t,i is considered to correctly model the pseudorange zero-mean measurement noise), m t,i = 0, r 2 t,i = 0.

(5.5)

• H 1 : the received measurement is in NLOS situation and affected by a mean value jump, m t,i = 0, r 2 t,i = 0.

(5.6)

• H 2 : the received measurement is in LOS situation and is affected by a variance change in the additive noise, m i,t = 0, r 2 i,t = 0.

(5.7)

Under hypothesis H 0 , the model error (i.e., the additive noise) has a Gaussian distribution. However, under hypotheses H 1 and H 2 , the nominal Gaussian distribution is no longer valid because of multipath presence. The errors associated to the two hypotheses H 1 and H 2 were modeled as Gaussian mixtures in [START_REF] Viandier | GNSS performance enhancement in urban environment based on pseudo-range error model[END_REF]. However, when analyzed more in detail, these mixture models (obtained from a real navigation scenario) can also be decomposed into a mean-shifted or a variance-increased Gaussian distribution (see Fig. 4 of [START_REF] Viandier | GNSS performance enhancement in urban environment based on pseudo-range error model[END_REF]). This decomposition considerably facilitates the analysis when compared to a Gaussian mixture model. Indeed, determining the number of components participating in a Gaussian mixture is fairly complicated.

System outline

Similarly to section 4.4.2 (eqs. (4.12)-(4.14)), the EKF innovations will be used to detect the presence of anomalies in the pseudorange measurements. A hierarchical method is proposed for error detection and later for error identification and correction. The idea is that this multi-stage approach enables an urban-adapted navigation filter without entailing heavy computations in clear sky scenarios. The mitigation scheme is described as follows:

1. Error detection:

The presence of an error is detected by performing a statistical test on the innovations

I t (H 0 ) = Y t -Y ( X t|t-1 ), (5.8)
where I t = (I t,1 , . . . , I t,ny ) and I t,i (H 0 ) ∼ N (0, s 2 t,i ) where s t,i is the ith diagonal element of the covariance matrix S t define in eq. (4.13) with r 2 t,i = 0.

Error identification:

In case an error has been detected, a parallel processing is achieved for classifying the two possible sources of error. Two "time of occurrence" tests are performed simultaneously for the detected outliers. These outliers can be affected by a mean value jump or by a variance change in the additive noise. The most likely hypothesis (H 1 or H 2 ) is then considered for error correction.

Error correction:

The innovation model is updated by correcting either the noise mean value or the noise variance, depending on the hypothesis that has been detected in the identification step. The corrected model is then fed back to the main system (composed by the EKF) that computes the final position.

Error detection

The first step of the algorithm detects the presence of corrupted signals referred to as outliers. The type of error affecting the signal is not specified at this point. This hierarchical approach, where the outliers are first detected and their source of error is later identified, is preferred to a strategy where the three hypothesis are simultaneously studied for every received signal because of its lower computational cost (considering that usually the number of received signals largely exceeds the number of outliers). A binary hypothesis test is performed to determine the absence (hypothesis H 0 ) or presence (hypotheses H 1 and H 2 ) of an error in the measurements. The test is achieved for each of the n y received signals through their respective innovations. A test based on the knowledge of the C/N 0 ratio was presented in [START_REF] Viandier | GNSS performance enhancement in urban environment based on pseudo-range error model[END_REF] to decide whether the received signal is error corrupted or not. However, if the multipath is in phase with the LOS signal, this test may no longer be valid. This chapter considers a sliding window of N samples as observation window and assumes that the error (when it exists) is constant during this period of time. The normalized energy of the innovations for each of the n y observation windows containing N samples is computed. The detection of errors is then achieved as follows:

T t,i = t j=t-N +1 I 2 j,i s 2 j,i H 0 ≶ H 1 or H 2 α ∀i = 1, . . . , n y (5.9)
where α is the detection threshold related to the probability of false alarm (PFA) of the test, and I j,i and s 2 j,i are obtained as explained in the previous section. The test statistics T t,i is distributed according to a central chi2 distribution with N degrees of freedom (denoted as χ 2 N ), under hypothesis H 0 . An accurate estimation of each innovation error distribution under H 0 (i.e. N (0, s 2 j,i )) is supposed to be available. The critical factor to obtain a precise value of s 2 j,i is to be able to correctly determine the nominal measurement noise variance σ 2 t,i associated to each ith measurement. Further details on the way this variance is calculated in practice are given in section 5.8.

If the test statistics T t,i exceeds the threshold α t,i , the presence of an error is declared for the ith measurement and the error estimation procedure is used to determine the kind of error affecting the received measurement. Fig. 5.4 shows hypothetic innovation error distributions under the hypotheses H 0 , H 1 and H 2 . It can be seen that the presence of error in the measurements (hypotheses H 1 and H 2 ) yields innovations with larger energy than under hypothesis H 0 . It is important to note that the threshold determination does not depend on the innovation distributions associated to H 1 and H 2 .

Moreover, no knowledge about the mean value jumps and the variance changes in the additive noise has to be known to compute the test statistics. The detected outliers are isolated from the solution in standard reliability tests. Such an approach is valid in clear sky scenarios where there is a large number of received measurements. However, in urban environments, visibility over the GPS constellation is scarce and measurement exclusion may lead to an undetermined system. Therefore the maximum number of received measurements is needed to compute the position solution. Thus the detection step must be followed by the estimation and correction of the anomaly. This is the objective of the following sections.

Error identification and correction

Once an outlier has been detected, its source of error has to be identified and corrected. Estimates of m and r are simultaneous computed and just the most likely type of error is used to determine the final navigation solution. Due to the recursive nature of the EKF, not only the error magnitude must be estimated but also its time of occurrence. Note that the error magnitude is generally smoothed within the observation window in (5.9), so that the effective time of detection of the error does not match its real time of occurrence k. In case of miss detection, the error gets propagated through the state vector, and the estimated error magnitude can differ significantly from the real value. As a consequence, both hypotheses H 1 and H 2 depend not only on the error parameter m or r (that are supposed to be constant inside the observation window) but also on the time of occurrence k (where k can take any value within the observation window containing samples from time instant t -N + 1

to time t) that should be estimated carefully. Change detection techniques such as the Neyman

Pearson test [VT01, page 33] need to know the model parameters conditionally to each hypothesis in order to provide an optimal solution. Unfortunately, this is not the case in practical applications where neither the time of occurrence of the change nor its magnitude are known. A test for error identification conveniently adapted to such situations is derived.

Time of occurrence estimation

The time of occurrence estimation can be achieved by using the marginalized likelihood ratio (MLR)

proposed by Gustafsson [START_REF] Gustafsson | Adaptive filtering and change detection[END_REF]. However, this test requires an a priori knowledge about the probability distributions of the parameters to be estimated (mean value jumps or variance changes in our case). An alternative is the generalized likelihood ratio (GLR) [START_REF] Willsky | A generalized likelihood ratio approach to the detection and estimation of jumps in linear systems[END_REF] which considers the jump as an unknown constant (as opposed to a random variable). The GLR proceeds for each qth detected outlier (where q ∈ [1, ..., n y ]) to a double maximization over the variable of interest (m or r) and the time of occurrence k. However, implementing the GLR in our system presents two main drawbacks.

On one hand, a bank of N no recursive least square filters (where n o denotes the number of detected outliers) must be deployed for each hypothesis (considering that the ensemble of outliers follows either H 1 or H 2 ), which presents a heavy computational cost for systems where the estimation delay must be extremely short. On the other hand, the GLR proposed in [START_REF] Willsky | A generalized likelihood ratio approach to the detection and estimation of jumps in linear systems[END_REF] is based on the idea that the relation between the error affecting the measurement and the EKF innovations can be made explicit and represented as a linear regression. However, even if a linear and recursive model can be easily stated for a mean value jump, this is not the case for noise variance changes. Hence, a time of occurrence test adapted to land navigation phenomena is proposed hereafter where just one filter is deployed for hypothesis. For the sake of simplicity, the estimation scheme will be developed under the assumption that the whole ensemble of n o outliers follows either H 1 or H 2 . It will be discussed later in this section the performance of the strategy in a more general situation where each outlier is independently affected by H 1 or H 2 .

A first estimation of the possible jump magnitudes ( m(k) or r(k), ∀k = t-N +1, . . . , t) is followed by the detection of the time of occurrence k. More precisely, the most likely time of occurrence k t,q

for hypothesis H 1 is defined as: k t,q (H 1 ) = arg min k {k|l t [k, m t,q (k)] > γ 1 } , (5.10) (5.11) where l t [k, m t,q (k)] is the log-likelihood ratio for the qth outlier under hypothesis H 1 at time instant k, and γ 1 is the test threshold that is fixed according to given prior probabilities for hypotheses H 1 and H 0 . Indeed, once the likelihood function overpasses the given threshold, k is chosen as the time of occurrence of the error (i.e. k t,q = k) and the associated mean value jump m t,q (k) is considered to remain constant for the rest of the observation window. The estimation of the mean value jump corresponding to the qth detected outlier and associated to a time of occurrence k is obtained as the mean value of the (t -k)th previous EKF innovation samples m t,q (k) = 1 t -k + 1 t j=k I j,q , ∀k = t -N + 1, ..., t.

(5.12)

where I j,q denotes the nominal qth innovation computed under hypothesis H 0 (i.e., I j,q = I j,q (H 0 )).

It should be noted that eq. (5.12) is valid under the assumption that a suitable observation window length has been chosen (i.e., the window is long enough to filter the noise contribution and short enough so that the state estimates have not significantly responded to the error). The log-likelihood ratio l t [k, m t,q (k)] for every possible time of occurrence k is calculated using the innovation samples at time instant k as follows l t [k, m t,q (k)] = log p (I k,q |H 1 ) p (I k,q |H 0 ) , (5.13) with p (I k,q |H 1 ) = p[I k,q (H 1 )], (5.14) where p(I k,q |H 1 ) represents the pdf associated to an error m t,q (k) at time instant k for I k,q , and where the EKF innovation I k,q (H 1 ) calculated under hypothesis H 1 is obtained from

I k (H 1 ) = Y k -Y [ X k|k-1 (H 1 )] -M k , (5.15)
where I k (H 1 ) is the innovation vector, I k,q (H 1 ) is its qth element, X k|k-1 (H 1 ) is the propagated state vector under hypothesis H 1 and M k is the mean value jump vector of dimension n y × 1. Each of the i elements of vector M k is defined as

• M k (i) = 0, if i = q (non corrupted measurement),
• M k (i) = m t,q (k), if i = q and the time of occurrence has not yet been detected for the qth outlier,

• M k (i) = m t,q ( k t,q (H 1 )), if i = q and the time of occurrence has already been detected for the qth outlier (i.e k > k t,q (H 1 )).

According to (5.15), each I k,q (H 1 ) follows a Gaussian distribution N (0, s 2 k,q ), where s 2 k,q is calculated under hypothesis H 0 as in (4.13). The propagated state vector X k|k-1 (H 1 ) in (5.15) is computed as X k|k-1 (H 1 ) = F X k-1 (H 1 ), where X k-1 (H 1 ) is the corrected state vector under hypothesis H 1 at time k -1. If no outlier has been detected at k -1 then X k-1 (H 1 ) = X k-1 (H 0 ). Hence, the corrected state vector at time k which depends on I k (H 1 ) and will be used at time k + 1 for the computation of X k+1|k (H 1 ), is defined as

X k (H 1 ) = X k|k-1 (H 1 ) + K k I k (H 1 ), (5.16)
where K k is the EKF gain matrix under nominal conditions (hypothesis H 0 ) and I k (H 1 ) is obtained from (5.15). Eqs. (5.15) and (5.16) highlight the straight relation between a correct innovation model and an unbiased position estimation contained in X k (H 1 ) and viceversa.

Similarly to the "time of occurrence estimation" for hypothesis H 1 presented in Eqs. (5.10) to

(5.14), the most likely time of occurrence under hypothesis H 2 is k t,q (H 2 ) = arg min k k|l t k, r 2 t,q (k) > γ 2 , (5.17) (5.18) where γ 2 is the test threshold that is fixed according to given prior probabilities for hypotheses H 2 and H 0 . The noise variance jump estimation r 2 t,q (k) for each possible time of occurrence k is defined as

r 2
t,q (k) = 1 t -k + 1 t j=k I 2 j,q -s 2 j,q .

(5.19)

The nominal innovation variance s 2 j,q (estimated under hypothesis H 0 and obtained from (4.13)) is subtracted from the calculated innovation variance in order to obtain the noise variance jump associated to the measurement. The likelihood function for H 2 is defined as

l t [k, r 2 t,q (k)] = log p (I k,q |H 2 ) p (I k,q |H 0 ) , (5.20)
for

p (I k,q |H 2 ) = p[I k,q (H 2 )], (5.21)
where p(I k,q |H 2 ) represents the pdf associated to an error r 2 t,q (k) at time instant k for I k,q , and where the EKF innovation I k,q (H 2 ) calculated under hypothesis H 2 follows a Gaussian distribution N (0, s 2 k,q (H 2 )) according to

R k (H 2 ) = diag( Θ k ), (5.22)
S k (H 2 ) = H k P k|k-1 H T k + R k (H 2 ), (5.23)
where s 2 k,q (H 2 ) corresponds to the qth element from the diagonal of S k (H 2 ) and Θ k is the measurement noise variance vector where each of its i elements is defined as

• Θ k (i) = σ 2 k,i , if i = q (non corrupted measurement),
• Θ k (i) = σ 2 k,i + r 2 t,q (k), if i = q and the time of occurrence has not yet been detected for the qth outlier,

• Θ k (i) = σ 2 k,i + r 2 t,q ( k t,q (H 2 )), if i = q and the time of occurrence has already been detected for the qth outlier (i.e k > k t,q (H 2 )).

The corrected state vector under hypothesis H 2 is then calculated as

X k (H 2 ) = X k|k-1 (H 2 ) + K k (H 2 )I k (H 2 ), (5.24)
where the EKF gain matrix K k (H 2 ) under hypothesis H 2 is computed using (5.23) according to

K k (H 2 ) = P k H T k (S k (H 2 )) -1 , (5.25)
and the innovation vector I k (H 2 ) (where I k,q (H 2 ) is its qth element) is described as follows

I k (H 2 ) = Y k -Y [ X k|k-1 (H 2 )].
(5.26)

In this way, Eqs. (5.22) to (5.26) highlight the dependence among the estimated noise variance parameter r 2 t,q , the final position solution contained in X k (H 2 ) and the innovations I k,q (H 2 ) to be used for the "time of occurrence test" in (5.20). Some pertinent remarks are made on the "time of occurrence test" algorithm related to both hypothesis H 1 and H 2 to facilitate the comprehension of the proposed approach:

1. For the computation of X k (•) all the available n y innovations are used. The proposed strategy was presented under the assumption that the whole ensemble of detected outliers was affected by the same type of error. For a more general case, where no correlation among the errors of the measurements is assumed, the straightforward strategy would be to build a bank of 2 no filters, each of them considering a possible combination of sources of error (mean value or noise variance jump) for the n o detected outliers. In urban areas where many simultaneous outliers can be present, this solution presents a high computational cost. However, we will now explain why the above developed strategy is also well suited to independent error measurements. Consider the general case where two measurements are simultaneously corrupted at time instant (t -d) (where d < N ) and that the error detection (5.9) is achieved at time instant t. Assume for instance that the first measurement Y 1 is affected by a mean value jump and the second measurement Y 2 is affected by a noise variance jump. At the output of each of the two filters considering exclusively H 1 or H 2 we will obtain:

• For the H 1 filter: once the H 1 test (5.10) is performed for the innovation samples corresponding to the (t -d) instant, the presence of an error will be detected in both measurements. For Y 1 a mean jump value will correctly compensate the innovation model. For the Y 2 innovations a mean jump compensation will also be achieved (because we are working under the H 1 filter).

However, according to (5.12), this estimated bias will approximately equal zero ( because Y 2 was just affected by a noise variance jump). As previously stated in the chapter, a non corrected noise variance jump affects mainly the covariance matrix of the state vector, not its mean value.

Therefore the Y 1 corrected innovations (5.15) won't be significantly affected by the non detected noise variance jump in Y 2 . This result is important for the next section, where the identification of the error is achieved through a likelihood test using all the compensated innovations of the observation window.

• For the H 2 filter: at (t -d) the H 2 test (5.17) will detect an error in the two measurements.

The innovation model corresponding to Y 2 will be correctly compensated by a jump in the noise variance. The Y 1 innovation model will be compensated by a fictitious noise variance jump.

The value of this variance jump will be proportional to the actual mean value jump affecting Y 1 (5.19). In theory, any non detected mean jump entails important biases in the state vector parameters (which will later condition the innovation calculation). However, the contribution of Y 1 to the state vector estimation is inversely proportional to its associated noise variance (i.e. if the noise variance is high, the EKF gain (5.25) for the measurement innovation will be low). Therefore, potentially high mean value jumps will be deweighted during the filtering.

In this way, the corrected Y 2 innovation model (5.26) won't be significantly affected by the non-compensated bias in Y 1 .

It is important to note that this discussion is valid even if the times of occurrence of the errors are not the same.

Error identification

After the time of occurrence detection was performed in parallel for H 1 and H 2 , a decision must be taken on the actual source of error affecting each outlier (i.e. either H 1 or H 2 ). Two different situations have to be considered at the output of the time of occurrence test:

1. An estimation of the error parameters corresponding to each outlier has been obtained under H 1 (with k t,q (H 1 ), m t,q ) and H 2 (with k t,q (H 2 ), r 2 t,q )

2. The presence of an error has been detected for only one hypothesis (i.e. for only one of the likelihood tests (5.10) or (5.17)).

In the second situation, no further test is required to decide on the measurement error source: the sole hypothesis under which the error was detected is considered as the actual source of error (mean value jump for H 1 or noise variance jump for H 2 ). However, if an error has been detected under both hypotheses, a final solution must be taken regarding the real origin of this error (remember that the measurement is exclusively affected by H 1 or H 2 ). Assuming both H 1 and H 2 have the same probabilities, the test for deciding between H 1 and H 2 (i.e. the measurements are affected by NLOS or LOS interference) is defined by p I t-N +1:t,q | k t,q (H 1 )

H 1 ≷ H 2
p I t-N +1:t,q | k t,q (H 2 ) (5.27) with p I t-N +1:t,q | k t,q (H 1 ) = t j=t-N +1 p [I j,q (H 1 )] ,

(5.28)

p I t-N +1:t,q | k t,q (H 2 ) = t j=t-N +1 p [I j,q (H 2 )] , (5.29)
where I t-N +1:t,q = (I t-N +1,q , ..., I t,q ).

Error correction

Once a decision has been taken regarding the source of error affecting each outlier, their associated innovation model has to be corrected and fedback to the update stage of the EKF algorithm. More precisely, if the qth outlier is affected by NLOS errors (H 1 ) the associated innovation is corrected according to

Ĩt,q = I t,q -m t,q k t,q (H 1 ) , (5.30) whereas if it is affected by LOS errors (H 2 ) the measurement variance is corrected as σ2 t,q = σ 2 t,q + r 2 t,q k t,q (H 2 ) , (5.31) where σ2 t,q and Ĩt,q denote the corrected parameters to be reinjected in the EKF algorithm. In the last case, the corrected measurement variance is used for the new computation of the measurement noise covariance matrix R t . In this way, a final unbiased navigation solution is calculated. It is important to observe that for the future error detection test in t + 1 (5.9), the Ĩt,q term will not replace the original innovation parameter I t,q calculated for H 0 . The observation window always contains the non corrected innovations to enable the detection of an error during its whole duration.

The final algorithm for the detection/identification/correction strategy is summarized in Table 5.1.

Results

Simulated Data

This section validates the proposed detection/estimation algorithm using simulated data. A land vehicle trajectory has been generated according to the state model (3.13) with acceleration standard deviations (stds) σ n = σ e = 2m/s and σ d = 0.2m/s (i.e. these acceleration stds apply to both the trajectory generator and the EKF's Q matrix). The received pseudorange measurements correspond to a simulated GPS constellation with nominal noise std σ = 12m. The number of visible satellites n y is constant and n y = 7. There are no changes in the constellation during the simulation period.

The satellite constellation is generated from online available GPS almanac data, with a 5 • visibility mask. The pseudoranges are supposed to be compensated by any type of atmospheric error. The errors introduced in the measurements have been generated according to the model (5.2) as follows • a mean value jump of 40m is introduced in the satellite number 1 between t = 30s and t = 60s,

• a noise std jump of 40m affects the same satellite for a time interval of 40 seconds between t = 100s and t = 140s,

• a second satellite (satellite number 2) experiences a mean value jump of 40m between t = 110s and t = 150s.

The simulated errors were introduced to highlight the performance of the proposed navigation filter.

The first isolated mean value jump on satellite 1 (corresponding to hypothesis H 1 ) is the type of error that more visibly affects the positioning accuracy. The correct functioning of the filter is tested for this critical situation. A simultaneous appearance of different errors is then studied. Two satellites are corrupted by different types of errors during overlapped time intervals. In this way, the algorithm is tested for its capacity to identify several defective measurements and their corresponding source of error. The threshold for the error detection in (5.9) has been adjusted in order to obtain P F A = 10 -5 . Based on several tests done using real data, a suitable observation window length of N = 5 was used. This choice was motivated by the need to detect fast changing errors and to achieve fast detection times, with no important losses in sensitivity. The data sampling period equals 1Hz. The estimation of the time of occurrence in (5.10) or (5.17) has been achieved with a threshold γ 1 = γ 2 = 1 (considering no a priori knowledge on the error time of occurrence, so both H 0 and H 1 or H 2 are given the same probability within the observation window).

The error detection/identification results for the interfered satellites are presented in Fig. The results for the final estimated position are presented in Fig. 5.9. The horizontal position errors (in 2D) are shown in blue and compared with their corresponding bounds illustrated in red.

The bounds are calculated as explained in section 4.7, for P c = 10 -5 . Result for a standard EKF (i.e.

Conclusions

This chapter presented an enhanced navigation system adapted to urban canyon scenarios. The originality of the proposed approach relies on the way the received signals are processed: a two step procedure is used to detect multiple outliers and to classify these outliers according to the different types of errors affecting the navigation signal. A hierarchical three-hypothesis test was implemented.

Two different situations were considered in the presence of multipath. These situations correspond to appeared as an interesting solution, due to the low quality of the exploited MEMS inertial sensors, results presented real weak performances. Differential odometry techniques are usually underestimated because their limitation to correctly compute the vehicle heading angle. However, results obtained within this thesis presented promising results for the GPS/WSS approach. A further strategy was tested using WSS and a 1-axis gyro especially developed for vehicle navigation. Though this gyro also belongs to the low-cost MEMS technology domain, its performances clearly outperformed those of the previous exploited inertial measurement unit. In this way, a final hybrid approach specially adapter to vehicle navigation could be proposed. Indeed, it is important to note that though an external gyro was used herein, most modern vehicles already contain an on-board gyro. Hence the WSS/Yaw-gyro DR strategy can be easily implemented at no additional cost. It must be also pointed out that the 1-axis gyro was not available at the beginning of the PhD, that is why its performances were only evaluated during the second stage of this work.

An alternative to the standard Extended Kalman Filter (EKF) used for navigation was studied.

The Unscented Kalman Filter(UKF) recently start gaining an increasing attention within the satellite navigation community. Its main advantage over the EKF is that no linearization must be applied to the non linear pseudorange measurement equation. However, in the particular context of urban land vehicle navigation using hybrid GPS/DR approaches the UKF didn't present any significant improvement over the standard EKF. This conclusion is in accordance with results presented in [START_REF] Gao | Development of precise GPS/INS/Wheel Speed Sensor/yaw rate sensor integrated vehicular positioning system[END_REF].

The presence of non-modeled measurement errors must be detected to avoid biased and inconsistent solutions. Typical errors affecting GPS measurement reliability are introduced by urban canyon phenomena in the form of multipath or cross correlation effects. Two different techniques were proposed with this purpose. High sensitivity strategies were applied in both cases to enable the tracking of typical low power urban signals.

GPS-based Doppler measurements are usually considered to be more robust to urban phenomena than pseudoranges. This is explained by the different techniques employed by the receiver to compute each of these measurements. Doppler measurements do not provide absolute position but only velocity estimations, and so they are not suitable for a stand-alone navigation solution. Consequently, a canyon phenomena present both on pseudorange and Doppler measurements.

• The choice of the confidence probability associated to the horizontal position bound computation should be standardized to account for any kind of urban phenomena. In this way, its implementation to safety of life applications in urban environments could be possible.

• It should be considered the application of the proposed multipath mitigation technique in chapter 5 to different areas where similar problems may be encountered. For example, similar errors were considered to affect the mobile communication signals [START_REF] Huerta | LOS-NLOS situation tracking for positioning systems[END_REF].

• The strategy proposed in chapter 5 could benefit from the use of Doppler measurements. We imagine an approach where Doppler information could be exploited to enhance the performance of the detection, identification and correction of erroneous pseudorange measurements. It will be also interesting to speculate about a fourth signal reception hypothesis where both a mean value and noise variance jump are considered.

xiv 4 . 9

 49 Errors for satellite 31. (a-b) pseudorange (PR) and Doppler measurement errors with their corresponding 3σ noise bound. (c-d) Histograms of normalized errors for pseudorange and Doppler measurements respectively. The nominal Gaussian PDF is represented in solid red line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Errors for satellite 5. (a-b) pseudorange (PR) and Doppler measurement errors with their corresponding 3σ noise bound. (c-d) Histogram of normalized errors for pseudorange and Doppler measurements respectively. . . . . . . . . . . . . . . . . . . . . . . . 4.11 Standard EKF-based navigation filter . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 RAIM+FDE strategy. Errors highlighted in green correspond to time instants where non reliable measurements were excluded while black dots indicate that RAIM+FDE strategy was not available. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.13 Navigation filter using the smooth correction criterion. Errors highlighted in green correspond to time instants where alternative measurements were used . . . . . . . . . 4.14 Number of alternative measurements used. Top figure: smoothed pseudorange (PR). Bottom figure: Doppler measurements . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.15 Results for reliability test using a binary correction criterion. In this case only Doppler measurements are exploited as an alternative to pseudoranges. . . . . . . . . . . . . . . 4.16 Vehicle circuit for different navigation systems in Toulouse centre. . . . . . . . . . . . . 4.17 Position errors (in blue) and bounds (in red) with GPS/WSS filter. (b) errors highlighted in green correspond to time instants where alternative measurements were used. 4.18 Position errors (in blue) and bounds (in red) with GPS/WSS/Yaw filter. (b) errors highlighted in green correspond to time instants where alternative measurements were used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Error envelope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Pseudorange errors in urban scenarios. (a-b) actual pseudorange errors for two different satellites. (c-d) normalized pseudorange error histogram (blue dotted line) and nominal Gaussian pdfs (solid red line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Proposed strategy for the detection, identification and correction of outliers. . . . . . . xv 5.4 Innovation pdf p(I t,i |H j ) (solid line for j = 0, dashed line for j = 1 and dashed-dotted line for j = 2). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Estimated pseudoranges obtained from the propagation of DR measurements. . . . . . 5.6 Error identification for the two interfered satellites. Dashed lines contain time intervals where errors are present. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.7 Innovation distributions for satellite 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Innovation distributions for satellite 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.9 Final position errors (in blue) and bounds (in red). . . . . . . . . . . . . . . . . . . . . 5.10 Position errors (in blue) and bounds (in red) with GPS measurements. For the error control strategies (b-c) and (d), the instants where a correction took place are shown in green. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Position errors (in blue) and bounds (in red) with GPS+DR measurements. For the error control strategy (b), the instants where a correction took place are shown in green.175 xvi List of Tables 1.1 GPS signal characteristics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Pseudorange error budget . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 The Extended Kalman Filter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 UKF algorithm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 UKF algorithm for additive noise. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 UKF scaling parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 HG1700 2D sensor noise parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 WSS sensor noise parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Asymptotic Horizontal PCRBs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 WSS simulated sensor noise parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Xsens MT9 2D sensor noise parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Horizontal error statistics in meters (error values highlighted in green represent the best performances). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Horizontal error statistics in meters for the GPS degraded intervals using the EKF . .

3. 12

 12 Horizontal error statistics in meters for GPS outage intervals . . . . . . . . . . . . . . 4.1 MLX90609 gyro noise parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Horizontal error statistics in meters and percentage of correct bounded error for GPS filter. Error values highlighted in green represent the best performances. . . . . . . . . 4.3 Horizontal error statistics in meters and percentage of correct bounded error for GPS/DR filter. Error values highlighted in green represent the best performances. . . . . . . . . xvii

1 .

 1 Augmentation of the GPS by using land vehicle DR sensors Micro Electro-Mechanical Systems (MEMS) technology has rendered inertial sensors suitable for civil land applications. Inertial navigation systems (INS) based on MEMS technology have become a standard research domain for land vehicle applications. However, due to relative lack of maturity of this technology, low cost INS show very limited performances. It is therefore interesting to study the upcoming of intelligent vehicle systems. The development of complex on board monitoring systems to control the vehicle trajectory and provide safety features entails the need for different DR sensors to be installed in the vehicle. In most new generation vehicles, options such as the Anti-lock breaking system (ABS) or the Electronic stability program (ESP) are becoming standard features. In this way, a new horizon in the vehicle land navigation domain is opened. Different sources of measurements describing the vehicle dynamics are available to be exploited at no additional cost. Wheel speed, tire pressure or the vehicle yaw angle are just some of information used by the typical advanced vehicle systems. The variety of available measurements presents an excellent opportunity to develop well-performing navigation systems.

Contents 1 . 1

 11 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 1.2 Satellite Navigation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.1 Global Positioning System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 1.2.2 Differential GPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.2.3 Urban navigation phenomena . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 1.2.4 High Sensitivity GPS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 1.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

  generating and distributing the navigation message to be uploaded into the satellite navigation message payload It is composed of three different sections aimed at performing the navigation tracking function, providing central command and control, and contacting and transmitting data to the satellites. Eleven monitoring stations (MS) are scattered around the globe and provide the Master Control Station (MCS) located in Colorado, with necessary raw pseudorange, carrier phase and metereological measurements. The MCS, which is the major brain of the control segment, processes this data to construct precise ephemerids and clock predictions. It also instructs the system if any change in the GPS satellite configuration must be achieved. The controls and updated navigation parameters are transmitted through the four exiting upload stations. The new navigation message is uploaded at least once a day to the satellite payload memory and transmitted to the GPS user.User segmentThis segment mainly deals with the GPS user receiving equipment. Any person equipped with a suitable GPS receiver around the world can process the received signals and obtained a position, velocity and time (PVT) solution. As already explained, the encrypted P(Y) code is only known by the military, so normal civilians won't have access to this high precision service. A commercial receiver will achieve the signal tracking for the signals with known C/A spreading codes.A block diagram for a generic digital receiver is presented in Fig.1.1. The received signal must first follow a "down-conversion stage" to allow all the posterior base-band processing. At this point the signal is ready to be processed by each of the N digital receiver channels where N determines the maximum number of possible satellites to be tracked (N is specific to each kind of receiver). Within these channels several functions such as loop discriminators and filters, data demodulation or C/N 0 computation are performed. Once the navigation message is demodulated, the ephemerids and error model parameters are processed to obtain the satellite position.

Figure 1 . 1 :

 11 Figure 1.1: GPS receiver block diagram

  t)) are the noises associated to the in-phase and quadrature signals respectively.The signal acquisition stage is required in order to look for a rough code delay and carrier frequency matching the received signal. This two dimensional search involves simultaneously finding a synchronized copy of the C/A code and carrier Doppler corresponding to the desired satellite. A search detector is deployed. For every delayed code/phase combination a statistical hypothesis test is computed. The test involves correlating the received signal with the local generated replica. If the two codes are aligned and the local carrier frequency matches the incoming signal's one, there is a maximum correlation. If the test exceeds a given threshold, the acquisition stage is considered completed.The tracking stage consists of performing a fine synchronization of the signal's phase and code delay considering the initial estimations obtained from the acquisition stage. A delay lock loop (DLL) is used to track the code delay while a phase lock loop (PLL) is used to track the carrier phase and the Doppler-shift due to the relative motion between the user and the satellite. Both tracking loops are simultaneously deployed and their output is used to generate synchronized local replica. If the alignment in the range domain (i.e. code delay domain) is successful, the code is removed from the signal so the navigation message can be demodulated. In parallel, the PLL generates sinusoidal signals to match the phase and frequency of the received signal. The carrier tracking loop may also be achieved by a frequency lock loop (FLL). Though this FLL is less sensitive to dynamic stress and therefore more robust to degraded scenario, it presents a lower accuracy. It is possible to used a FLL-assisted-PLL carrier tracking loop to improve performances.

GPS navigation in urban

  scenarios (contrary to open-sky scenarios) presents one of the most challenging topics within the GNSS community. The presence of several urban-based error sources such as multipath, interference, signal masking and poor constellation geometry entails an important loss of precision in the navigation solution. Signal attenuation is caused by the presence of obstacles such as trees in the LOS path between the user and the receiver. Weak signals arriving to the receiver seriously degrade the navigation performance. Cross-correlation errors appear in the presence of noisy CHAPTER 1. GPS SATELLITE NAVIGATION signals.

Figure 1 . 2 :

 12 Figure 1.2: Correlation function in the presence of one multipath component. The final correlation function in full red line is the result of the addition between the direct signal (dashed blue line) and multipath (dash-dotted green line) correlation function.

Figure 1 . 3 :

 13 Figure 1.3: Modes of multipath presence in urban environments.

Contents 2 . 1

 21 Dead reckoningDead reckoning (DR) is the use of sensors that provide relative positioning information[START_REF] Bullock | A prototype portable navigation system utilizing map aided GPS[END_REF]. DR navigation is based on the knowledge of the vehicle dynamics and an original starting point to calculate the actual vehicle position. Different sensors provide attitude rate angle and velocity measurements.These measurements are used to propagate the position estimations done on previous time instants and thus keep track of the vehicle trajectory.

Figure 2 . 1 :

 21 Figure 2.1: Accelerometers scheme [Ste00].

Figure 2 . 2 :

 22 Figure 2.2: Idle sensor output [Cal07].

Figure 2 . 3 :

 23 Figure 2.3: Typical Allan deviation plot for a system [Cal07].

Fig. 2 .Figure 2 . 4 :

 224 Fig. 2.4 depicts a WSS based on the hall-effect [Ste00]. A ferrous wheel (i.e. the sensor ring) is attached to the vehicle's wheel. When a current is passed through the hall-effect device present in the so called WSS, a voltage develops across the device perpendicularly to the direction of the current flow and to the direction of magnetic flux. As the sensor ring rotates, the reluctance of the magnetic field changes as the teeth pass the sensor. By multiplying the wheel rotational rate (in rad/s) by the

Figure 2 . 5 :

 25 Figure 2.5: Speed effect on tire diameter [Ste00].

Figure 2 .

 2 Figure 2.6: Device to recover the WSS data.

Figure 2

 2 Figure 2.7: Main coordinate frames used in inertial navigation.

  herein, the mobile axis m = (X m , Y m , Z m ) are aligned with the inertial measurement unit (IMU) frame. Its orientation with respect to the navigation frame is defined by the Euler angles (φ, θ, ψ) standing for roll, pitch and heading (or yaw) angle. Each of these angles results from a rotation around one of the mobile axis. In particular, as a 2D solution will be considered in this thesis, only the rotation around the Z m axis (i.e. the heading angle) is of interest. Fig.2.8 depicts the mobile frame considering just a heading rotation.

Figure 2 . 8 :

 28 Figure 2.8: Mobile frame (m).

2 .

 2 Transformation of specific force to navigation frame of interest 3. Velocity calculation4. Position calculationThe attitude update is obtained from (2.21) and the frame transformation of the accelerometers' specific forces is given by the matrix R m2n presented in (2.16). Velocity and position calculations are respectively given by (2.15) and (2.12). A schematic illustration of the 2D INS mechanization, neglecting centrifugal forces and considering the approximation given in (2.20), is depicted in Fig.2.9.

Figure 2 .

 2 Figure 2.9: INS mechanization.

  rameters are considered. The final INS position solution is given by (λ IN S + δλ, ϕ IN S + δϕ), where (λ IN S , ϕ IN S ) is the calculated INS position (from the INS mechanization) without any error correction.

Figure 2 .

 2 Figure 2.10: Illustration of WSS definition.

  and velocity estimation. The IMU outputs are fed to the INS computer to get position, velocity and attitude estimation. GPS and INS navigation parameter estimations are then fed to the "coupling or integration filter" as shown in Fig. 3.1. The updated state covariance matrix of the GPS filter is commonly used as the measurement covariance matrix in the coupling filter. The DR solution is finally corrected by the estimated DR errors.

Figure 3 . 1 :

 31 Figure 3.1: Loose coupling integration scheme.

Figure 3 . 2 :

 32 Figure 3.2: Tight coupling integration scheme.

Figure 3 . 4 :

 34 Figure 3.4: Unscented transform.

  .33) where i = 1, ..., n x and n x is the dimension of the state vector. The parameter k is a scaling parameter and (n x + k)P i denotes the ith column (or row) of the matrix square root of the weight covariance matrix (n x + k)P . The matrix square-root is usually obtained by means of the Cholesky factorization. The values for the weights and sigma point locations are obtained by comparing the Taylor series expansions of the estimators (3.27)-(3.29) with the true quantities and choosing the weights and sigma

Fig. 3 .

 3 5(b) shows the positioning errors for the GPS/WSS coupling approach which are compared to the corresponding PCRB.These errors fluctuate from 10 to 15 meters. Note that the coupling technique performs similarly when using the EKF or the UKF. The vehicle dynamics is sufficiently slow for this example to linearize efficiently the measurement equation (3.52). Therefore the UKF doesn't show any improvement over the EKF. Note also that the hybridization approach reaches almost the best achievable precision, provided by the Posterior Cramer Rao bound. Similar results obtained for the GPS/INS/WSS hybridization can be observed in Figs. 3.6(a) and 3.6(b). A comparison between both coupling techniques shows that the GPS/INS/WSS outperforms the GPS/WSS approach. Adding INS measurements to GPS and WSS data results in improved accuracy. For instance, the GPS/INS/WSS positioning errors fluctuate from 5 to 10 meters.

Fig. 3 .

 3 Fig. 3.7 compares the positioning errors obtained for the different coupling and filtering strategies for this example. The GPS/INS/WSS coupling clearly outperforms GPS/WSS during limited GPS visibility and GPS outages. The use of INS helps to keep a bound over the positioning errors, avoiding loosing track of the vehicle. Consequently, GPS/INS/WSS hybridization provides reliable positions during longer periods than the GPS/WSS system. Note also that a slight improvement is observed with the UKF (versus the EKF) when regaining full GPS visibility. Another way of comparing the different coupling strategies is to analyze the corresponding horizontal position PCRBs for different numbers of LOS satellites. Table3.7 shows that the asymptotic

Figure 3 . 5 :

 35 Figure 3.5: GPS/WSS approach.

Figure 3 .Figure 3 . 7 :

 337 Figure 3.6: GPS/INS/WSS approach.

Figure 3 . 8 :

 38 Figure 3.8: Vehicle circuit in ISAE campus, Toulouse.

3 .

 3 11 for the yaw-rate measurement. Figs. 3.10(a) and 3.11(a) depict the resulting error sequences and Figs. 3.10(b) and 3.11(b) illustrate the respective error histograms. The chosen noise variances 100CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA Error histogram for WSS velocity

Figure 3 .

 3 Figure 3.10: Characterization of the noise variance for the WSS velocity.

  Figure 3.11: Characterization of the noise variance for the WSS yaw-rate.

102CHAPTER 3 .

 3 HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

Figure 3 .

 3 Figure 3.12: Allan variance for the Xsens MT9 gyroscope.

Figure 3 .

 3 Figure 3.13: Allan variance for the Xsens MT9 accelerometers.

Figure 3 .

 3 Figure 3.14: Horizontal error (HE) for stand alone GPS.Performances for the different hybrid approaches are presented in Fig.3.17 and summarized in Table3.10 and Table3.11. To account for error statistical properties (which are not as straightforward as for the simulated data), different parameters are computed. The 50% bound represents a lower threshold under which half the error values are found. The bound is obtained by arranging the errors

Figure 3 .

 3 Figure 3.15: Number of visible satellites.

Figure 3 .

 3 Figure 3.16: Intervals subjected to measurement errors.

Figure 3 .

 3 Figure 3.17: Horizontal error for different filters and hybrid systems.

Fig. 3 .

 3 Fig. 3.18(a) presents the location of the outages in the vehicle trajectory. The locations were chosen to test the DR performances for different dynamics (straight or curved trajectory). The respective number of visible satellites is depicted in (Fig. 3.18(b)). Satellites with good elevation angles (>60 • )

Figure 3 .Figure 3 .

 33 Figure 3.18: GPS outage characterization.

Figure 3 .

 3 Figure 3.20: Stand alone WSS navigation solution.

Figure 3 .

 3 Figure 3.21: Stand alone INS solution.

  the least square case, and of the innovation sequence in the case of the Kalman filter. The detection of model errors is based on statistical hypothesis testing. Assuming a correct nominal measurement model, observation residuals/innovations indicate the extent to which the measurements are internally consistent with their assumed error model. Usual strategies employed by commercial receivers to verify the quality of the received signals are based on the receiver autonomous integrity monitoring (RAIM) strategy [C + 06a]. This strategy is based on the least squares observation residuals and it is presented in section 4.4.1. It is usually complemented by a fault detection estimation (FDE) technique. The RAIM+FDE strategy is presented herein as a test-bench solution. Section 4.4.2 presents an innovative strategy for mitigating the presence of unreliable pseudorange measurements. Doppler measurements are proposed herein as an alternative to fight against urban phenomena.

  ζ the nominal hypothesis is kept (in our case the pseudorange measurement is exploited) and if T t,i > ζ an anomaly in the measurement is declared (in our case the pseudorange is replaced by a Doppler measurement). The correct choice of the threshold is one of the most delicate issues for the well-conditioning of these kinds of tests. Indeed, we propose to relax the critical role of the threshold computation by replacing it by an interval. Two thresholds (ζ 1 , ζ 2 ) are now related to the test statistics. An extra degree of freedom is thus introduced to the global anomaly detection. Test results performed for each of the ith visible satellites can lead to three different conclusions:

• ζ 1

 1 < T t,i < ζ 2 the smoothed pseudorange measurement presented in (4.7) is used. The merging weight α t,i is determined by a smooth correction function which linearly depends on the values taken by the test between the interval ζ 1 < T t,i < ζ 2 .

Figure 4

 4

Figure 4 . 1 :

 41 Figure 4.1: Test statistic for the EKF innovation.

Figure 4 . 2 :

 42 Figure 4.2: Weight function.

Figure 4 .

 4 Figure 4.3 summarizes the above mentioned constraints applied to the ith received signal. Considering that urban canyon phenomena are correctly mitigated by the use of smoothed pseudorange and Doppler measurements, the final position estimation should be consistent and accurately described by the updated EKF state covariance matrix (table 3.1).

Figure 4 . 3 :

 43 Figure 4.3: Constraints on the choice of the type of measurement to be exploited.

Figure 4 . 4 :

 44 Figure 4.4: Horizontal position error ellipse.

A

  deep urban trajectory was chosen to analyze the relevance of the proposed smooth correction criterion. It is important to notice that this chapter fully focuses on the challenges and impairments of processing real data. Assumptions made on the behavior of errors affecting the received signals need to be analyzed as well as pseudorange and Doppler characteristics. Therefore, using simulated data to obtain theoretical performances would go against the intrinsic "uncertain and random" behavior proper to urban phenomena. A test field campaign was done in Toulouse centre for a duration of 25 minutes. The total distance traveled was 5.6km. Fig.4.5 depicts the vehicle trajectory and Fig.4.6 shows the number of visible satellites as a function of time. The scenario was carefully chosen to highlight the impact of urban effects. Specifications for the reference position solution and used GPS receiver were given in section 3.8.2. ABS WSSs are characterized according to noise models given in section 3.8.2. The 1 axis MLX90609 gyro proposed for the GPS/WSS/Yaw approach is characterized in section 4.8.1. An analysis on the theoretical robustness of Doppler measurements to urban phenomena is done in section 4.8.2. Results on navigation performances are given in section 4.8.3.

Figure 4 . 5 :

 45 Figure 4.5: Vehicle circuit in Toulouse centre.

Figure 4 . 6 :

 46 Figure 4.6: Number of visible satellites.

Figure 4 . 7 :

 47 Figure 4.7: Allan variance for the Melexis gyroscope.

Fig. 4 .

 4 8 depicts the satellite constellation. In particular, a good elevation (satellite 31) and a low elevation (satellite 5) satellite are chosen to illustrate a general case. The SPAN reference solution is used to estimate the errors affecting the received measurements.Theoretical pseudorange and Doppler measurements are reconstructed according to their model given in (3.50) and (4.2) excluding the noise term. The vehicle's position and velocity are provided by the reference solution while the clock bias and drift are obtained from the GPS Ublox receiver. In this way, any difference between the theoretical measurements and the real ones will account for noise and urban phenomena.

Figs. 4 .

 4 Figs.4.9 and 4.10 present pseudorange and Doppler measurement errors for satellites 31 and 5

Figure 4 . 8 :

 48 Figure 4.8: Satellite constellation.

Figure 4 . 9 :

 49 Figure 4.9: Errors for satellite 31. (a-b) pseudorange (PR) and Doppler measurement errors with their corresponding 3σ noise bound. (c-d) Histograms of normalized errors for pseudorange and Doppler measurements respectively. The nominal Gaussian PDF is represented in solid red line.

Figure 4 .

 4 Figure 4.10: Errors for satellite 5. (a-b) pseudorange (PR) and Doppler measurement errors with their corresponding 3σ noise bound. (c-d) Histogram of normalized errors for pseudorange and Doppler measurements respectively.

Figure 4 .

 4 Figure 4.13: Navigation filter using the smooth correction criterion. Errors highlighted in green correspond to time instants where alternative measurements were used

Figure 4 .

 4 Figure 4.14: Number of alternative measurements used. Top figure: smoothed pseudorange (PR).

4 .

 4 16 illustrates the stand-alone DR positioning solution for the WSSs and WSS/Yaw strategies. The accurate gyro clearly helps to overcome WSSs limitation in the heading angle computation.

Figure 4 .

 4 Figure 4.16: Vehicle circuit for different navigation systems in Toulouse centre.

Figs

  Figs.4.17 and 4.18 present the advantages of implementing the smooth correction criterion. While

  et al.[START_REF] Giremus | Joint detection/estimation of multipath effects for the global positioning system[END_REF] studied a Rao Blackwellized particle filter based on a jump Markov system. The proposed algorithm modeled the multipath NLOS situation by a mean value jump whose magnitude was jointly estimated with the vehicle position and velocity. Another two-hypothesis Bayesian approach was considered in[START_REF] Viandier | GNSS performance enhancement in urban environment based on pseudo-range error model[END_REF]. The interfered signals were characterized by error models based on

Figure 5 . 1 :

 51 Figure 5.1: Error envelope

  5.2 presents results for the test field campaign described in section 4.8. Similarly to section 4.8.2, the pseudorange errors for two different satellites are shown in the top figures. The histograms of the normalized errors (i.e. errors divided by their corresponding noise standard deviation) are displayed in the bottom figures. In open sky scenarios the normalized pseudorange errors follow a Gaussian probability density function (pdf) [Lei04]. Conversely, the bottom figures show that the Gaussian pdf assumption is not valid for the considered satellites. In Fig. 5.2(c) the noise variance seems to be underestimated by the fitted Gaussian pdf, whereas Fig. 5.2(d) shows a mean shifted error histogram. Note that the pseudorange errors presented similar characteristics in [VNMD08] and [DH08] under different conditions.

Figure 5 . 2 :

 52 Figure 5.2: Pseudorange errors in urban scenarios. (a-b) actual pseudorange errors for two different satellites. (c-d) normalized pseudorange error histogram (blue dotted line) and nominal Gaussian pdfs (solid red line).

Fig. 5 .

 5 Fig. 5.3 illustrates the above mentioned approach. Details about the detection, estimation and correction of errors are provided in the following sections.

Figure 5 . 3 :

 53 Figure 5.3: Proposed strategy for the detection, identification and correction of outliers.

Figure 5 . 4 :

 54 Figure 5.4: Innovation pdf p(I t,i |H j ) (solid line for j = 0, dashed line for j = 1 and dashed-dotted line for j = 2).

2 .

 2 Only the innovations associated to the detected n o outliers are corrected.3. The innovation vector Ik (•) contains (n y -n o ) non corrected innovations.

Figure 5 . 5 :

 55 Figure 5.5: Estimated pseudoranges obtained from the propagation of DR measurements.

  5.6. The top figure shows results for satellite 1 and the bottom figure for satellite 2. The strategy reveals a very good performance where the correct error hypothesis (H 0 , H 1 or H 2 ) is almost always detected.

Figs. 5 .Figure 5 . 6 :Figure 5 . 7 :

 55657 Figs. 5.7 and 5.8 show the innovation pdfs corresponding to the two corrupted satellites. The nominal Gaussian pdf is depicted in red while the actual normalized pdf for the EKF innovations is shown in blue. The pdfs are obtained from all the available samples of the simulated satellites. Top figures (a) present results for a standard EKF while bottom figures (b) correspond to the enhanced algorithm proposed in this paper. The innovations do not have a Gaussian distribution in the first case because the corrupted measurements are not compensated by the filter. Conversely, when errors are corrected with the proposed algorithm, the histograms of the corrected innovations are close to the adjusted Gaussian pdf.

Figure 5 .

 5 Figure 5.11: Position errors (in blue) and bounds (in red) with GPS+DR measurements. For the error control strategy (b), the instants where a correction took place are shown in green.

Table 1 .

 1 1 summarize the GPS signal characteristics.

	GPS band	L1		L2
	Carrier Frequency	f 1 = 1575.42M Hz	f 2 = 1227.60M Hz
	Code	C/A	P(Y)	P(Y)
	Power	-160dBW	-163dBW	-166dBW
	Chip rate	1.023MHz 10.23MHz	10.23MHz
	Period	1ms	1 week	267 days
	Navigation message	Data rate Chip width 50 bps 20 ms	Period 12.5 min

Table 1 .

 1

1: GPS signal characteristics C/A codes are widely known and available to any civilian user, while the P(Y) codes are only available to the U.S military and other authorized users. This limited access to code information entails to different positioning services: the civil Standard Positioning Service (SPS) and the military Precise Positioning Service (PPS). The service provided on L2 and L1 with the P(Y) code is not considered within this thesis. All the results present hereafter will correspond to measurements obtained from the SPS. The PRN codes are used for spreading the signal frequency spectrum and to identify the received satellite. This CDMA (code division multiple access) technique allows several signals to share the same spectrum bandwidth with limited interference among each other. The original signal power is spread through a large bandwidth so that it gets almost masked by thermal

  Filtering framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 The optimal Kalman Filter (KF)-Linear systems . . . . . . . . . . . . . . 3.4 Non-linear systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Cramer Rao bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 GPS/Dead-reckoning-sensors hybrid systems . . . . . . . . . . . . . . . . . Discrete state model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.1 Synthetic data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.8.2 Tests done on real data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

GPS, INS and odometric data

Contents

3.1 Integration strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Loose integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.2 Tight integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Chosen integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 3.4.1 The Extended Kalman Filter (EKF) . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 The Unscented Kalman Filter (UKF) . . . . . . . . . . . . . . . . . . . . . . 3.5 3.6.1 GPS/WSS integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6.2 GPS/INS/WSS integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7

•

  Update the state mean and covariance matrix with the received observations Xt = Xt|t-1 + K t I

t P t = (I nx -K t H t ) P t|t-1

Table 3 .

 3 1: The Extended Kalman Filter.

Table 3 .

 3 3: UKF algorithm for additive noise. 80CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

  INS error state model describing F IN SOriginally, the GPS/INS systems were aimed at aerial and spatial applications where highly precise IMU were used. However, for commercial land vehicle navigation such inertial systems are not affordable and thus lower quality sensors must be used. The problem is that during GPS outages the errors of the INS for this type of low cost IMU can reach unacceptable values. For a low cost IMU as the one employed in this thesis with sensors based on micro electro-mechanical system (MEMS) technology, a positioning error of 100 meters can be attained after 30 seconds of GPS outage[START_REF] Godha | Performance evaluation of low cost MEMS-based IMU integrated with GPS for land vehicle navigation application[END_REF].

	86CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA
	The measurement model will be analogous to (3.52) where the position errors are now obtained from
	the state vector in eq. (3.56) and are used to correct the INS position solution as (λ IN S t	+ δλ t , ϕ IN S t	+
	δϕ t , h t ).
									t	T GPS/INS/WSS	∈ R 12 .	(3.58)
	Considering (3.57), the linear state transition matrix	F GP S/IN S/W SS t	and the state noise matrix
	BGP S/IN S/W SS
					t	and BINS t	was already given at the end of chapter 2, so
	no further discussion will be done herein. Similarly to (3.48), the state matrices are constructed as:
	F GP S/IN S t	=	 	F IN S t 0	0 F GP S t	  ,	BGP S/IN S t	=	 	BINS t 0		0 t BGP S	  .	(3.57)

To overcome this loss of accuracy, a second source of dead reckoning measurements is going to be introduced. In particular we propose to exploit the already on board odometric data provided by the ABS wheel speed sensors.

GPS/INS/WSS

State model

The state vector (3.56) is augmented by radius errors affecting the WSS measurements yielding Ẋt = δ ṗt , δ vt , δ ψt , ḃa,t , ḃg,t , ḃt , ḋt , δ Ṙrr t , δ Ṙrl t for the GPS/INS/WSS approach in a time-continuous space are obtained as:

  The expression for h ψ (.) is given by the difference between the INS angle error δψ and the radius errors propagated through the WSS yaw

	90CHAPTER 3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA
	and the WSS-based measurements ((3.60) and (3.70)) as:
	δR rr t , δR rl t , δψ t ) + n ψ t ,		(3.70)
	where (ψ rate expression (3.45):
	h ψ (δR rr t , δR rl t , δψ t ) = δψ t -	ω rr t δR rr t -ω rl t δR rl t L	.	(3.71)
	The final measurement vector is thus given by the concatenation of the GPS measurements (3.49)

t ) W SS is the yaw angle measured by the WSSs, (ψ t ) IN S is the yaw angle measured by the INS and n ψ t ∈ R 1 is the noise vector n ψ t ∼ N (0, σ 2 ψ ).

  3. HYBRID NAVIGATION SYSTEMS: INTEGRATING GPS, INS AND ODOMETRIC DATA

			Mean 50% bound 95% bound
		GPS	6.85	6.63	13.27
	EKF	GPS/WSS	7.44	6.71	13.51
		GPS/INS/WSS 6.74	6.67	13.15
	UKF	GPS/WSS GPS/INS/WSS 6.78 7.41	6.73 6.73	13.68 12.81
	Table 3.10: Horizontal error statistics in meters (error values highlighted in green represent the best
	performances).
		1 st interval with degraded pseudoranges
			Mean 50% bound 95% bound
		GPS	12.20	9.66	21.95
		GPS/WSS	9.69	8.70	13.32
	GPS/INS/WSS 8.26	6.95	13.41
		2 nd interval with degraded pseudoranges
		GPS	13.03	13.63	18.12
		GPS/WSS	9.65	8.62	15.23
	GPS/INS/WSS 9.44	8.56	17.51
	Table

  3.8. RESULTSsimilar performances. For both partial or total outages the WSS-based strategy clearly outperforms the INS/WSS-based strategy. In contrast with the simulated scenario, where a better quality IMU was simulated, the Xsens MT9 IMU doesn't show any positive contribution to the navigation solution.In the presence of partial or total outages, the available measurement (pseudoranges, if present, and WSS measurements) are not enough to estimate all the errors of the high dimension INS state vector. Hence, no serious control can be obtained over the exponentially growing INS position errors.Either a biased (partial outage) or diverging solution (total outage) is obtained. However, WSSbased navigation solution is still reliable for GPS outages. Even for the 2 nd total GPS outage, which takes place during a 360 • turn, the WSS errors remain bounded. Fig 3.20 and Fig. 3.21 emphasize the contrast in the performances obtained by the stand-alone WSS-based navigation solution when compared to the stand-alone INS. It is clear that no substantial enhancement can be introduced by

		0.7604
		0.7604
	Latitude [rad]	0.7604 0.7604 0.7604	15sec, partial outage				20 sec, total outage
				20 sec, total outage
		0.7604
		0.7604
		0.7604
			0.0257	0.0257	0.0257	0.0257	0.0257	0.0257	0.0257
					Longitude [rad]

the low-cost Xsens MT9 IMU , if WSS data allowing for differential odometry (i.e. velocity and yaw-rate measurements) is available.

  Ẋt = λt , φt , vn,t , ve,t , ḃt , ḋt given in the previous chapter. The state propagation matrix F t corresponds to a random walk model where the acceleration has a zero mean Gaussian distribution modeled by the state noise covariance matrix Q t . The acceleration variances associated to (v n,t , v e,t ) are denoted as σ 2 n , σ 2 e . In a continuous time model the state matrices are constructed as follows:

	4.5 State-space models for GPS filter using smooth correction crite-
	rion
	In a first approach, a stand-alone GPS filter is used to test the relevance of the proposed smooth
	correction criterion. The time continuous state vector (always assuming a 2D navigation with known
	height) is given by:
	T GPS	∈ R 6 ,	(4.18)
	according to notations

  The filter state vector X Ẋt = δ λt , δ φt , δ Vt , δ ψt , δ Ṙrl

	GPS/WSS/Yaw t	, where Yaw denotes the yaw gyroscope, is given
	by,
		T GPS/WSS/Yaw	∈ R 9 .	(4.26)
	where
	• (δλ t , δϕ t ) are the positioning errors corresponding to the DR navigation solution (λ WSS/Yaw t	WSS/Yaw , ϕ t

t , δ Ṙrr t , ḃg,t , ḃt , ḋt

  WSS/Yaw

			t	=	 	BWSS/Yaw t 0	0 t BGPS	  ,	(4.28)
	where	BWSS/Yaw t	differs from the BWSS t	matrix given in (2.31) in that the yaw angle noise variance
	σ 2 ψ is now determined by the gyro bias variance σ 2 bg . The expression for BGPS

t was defined in

(3.46

).

Table 4 .

 4 1: MLX90609 gyro noise parameters.

Table 4 .

 4 2: Horizontal error statistics in meters and percentage of correct bounded error for GPS filter.

		Mean 50% bound 95% bound Bounded error
	Standard GPS filter	8.77	6.97	19.88	90%
	Smooth correction criterion 7.82	6.38	19.55	96%
	Binary correction criterion	8.62	7.67	17.7	63%

Table 4 .

 4 . 4.17 and 4.18 present the advantages of implementing the smooth correction criterion. While 3: Horizontal error statistics in meters and percentage of correct bounded error for GPS/DR filter. Error values highlighted in green represent the best performances.

	Mean 50% bound 95% bound Bounded error

  ) where • Y t,i is the ith pseudorange measurement associated to the ith visible GPS satellite for i = 1, ..., n y (n y being the number of visible satellites), ) is the ith GPS satellite position expressed in ECEF rectangular coordinates,• p t is the vehicle position expressed in ECEF rectangular coordinates. It is obtained from the state variables (λ t , φ t , h

	• p s t,i = (x s t,i , y s t,i , z s t,i

t , considering h t known) with an appropriate frame transformation (for more details see [FB99, ch. 2]),

Table 5 .

 5 2: Horizontal error statistics in meters and percentage of correct bounded error. Error values highlighted in green represent the best performances.

		Mean 50% bound 95% bound Bounded error
	GPS	Standard DIC approach 7.51 8.77	6.97 5.82	19.98 18.01	90% 98%
	GPS/WSS/Yaw	Standard DIC approach 6.18 7.73	6.29 5.53	18.52 12.77	91% 96%

Acknowledgements

Chapter 4

Doppler smoothing

Contents CHAPTER 4. DOPPLER SMOOTHING solutions were successfully corrected. This means that the pseudorange outliers were correctly detected and replaced by suitable smoothed pseudorange or Doppler measurements. The best performances in terms of mean and bounded error were obtained for this strategy.

• Fig. 4.14 outlines the contribution of the smoothed pseudorange measurement. The 58% of the times a pseudorange outlier was detected the smoothed pseudorange was used, against the 42% of the times for the Doppler measurement. Hence, the non corrected peaks in Fig. 4.13. Therefore, if Doppler measurements are not correctly exploited, they can produce biased solutions. Consider this situation takes place at time instant t.In the following time instant t + 1 a false outlier detection will be declared preventing pseudoranges from compensating the position error. Even if in t + 1 the Doppler measurement is no longer erroneous, it will only provide velocity information to propagate the previous biased position estimation. Hence, the contribution introduced by the smooth correction criterion where an extra degree of freedom is introduced to relax the fault detection threshold computation. parameter is aimed at providing information on the percentage of horizontal position errors that are CHAPTER 4. DOPPLER SMOOTHING the use of DR sensors already helps to filter the presence of "measurement errors" (Figs. 4.17(a) and 4.18(a)), the inclusion of the smooth correction criterion enables very good performances to be obtained (Figs. 4.17(b) and 4.18(b)). Results for the GPS/WSS and GPS/WSS/Yaw approach are summarized in table 4.3.

Multipath mitigation by error identification

Contents

Initialization

• PFA (used for the calculation of α t,i and γ), window length N and state matrices X 0 and P 0 .

For t = 1, . . . , end:

• Initialize: Ĩt,i = I, i t and σ2 t,i = σ 2 t,i , ∀i = 1, . . . , n y .

• If t > N (enough samples in the observation window):

• Detection: error detection test according to (5.9)

Implementation of an EKF considering hypothesis

Computation, for each outlier, of the likelihood ratios in (5.13) using m t,q (k) until the time of occurrences k t,q (H 1 ) are detected according to (5.10).

Computation of the innovations I k (H 1 ).

Implementation of an EKF considering hypothesis H 2 for k = t -N + 1, . . . , t.

Computation, for each outlier, of the likelihood ratios in (5.20) using r 2 t,q (k) until the time of occurrences k t,q (H 2 ) are detected according to (5.17).

Computation of innovations I k (H 2 ) and their associated variances s 2 t,q (H 2 ).

Identification: compute likelihood tests in (5.27) to identify H 1 or H 2 .

Correction: the innovation model parameters for each outlier are corrected according to their identified source of error, replace

, and

Table 5.1: Detection/identification/correction strategy to mitigate multipath interference.

Land navigation system with GPS/WSS/Yaw hybrid approach

State-space models for the GPS/WSS/Yaw approach were presented in section 4.6 (where just the pseudorange measurement model is considered herein). Though the reliability test principle (as already detailed in sections 5.5 and 5.6) does not change for the augmented system, a slight modification is introduced in the computation of the observation window innovations I j,i used in (5.9), (5.12) and

(5.19). The idea is to exploit the immunity of the DR (i.e. WSS/Yaw) navigation system against multipath effects to compute these innovations. The N innovation samples corresponding to the observation window are not taken from the already calculated innovations in (5.8), where the state vector was used to predict the received measurements. Instead, they are recalculated using only the DR measurements as follows

where obtained under H 0 , but is iteratively obtained as,

with

where β is initialized as β t-N = P t-N , P t-N being the updated state covariance matrix of the EKF at time instant t -N . The estimated pseudorange measurements Ỹj,i are obtained according to (5.33) by propagating an initial position estimation a t-N with ∆a DR u . Therefore, a faster and more accurate error detection is achieved. However, it must be observed that this DR position propagation strategy is independent of the EKF implementation, i.e., the new I DR j,i does not replace the standard EKF innovations (5.8) at time j and is only used for (5.9), (5.12) and (5.19). This approach presents the advantage of including a second navigation system, not affected by multipath, that enables a more efficient detection/identification/correction strategy. In the following section, the performance of the augmented DR system is tested with real navigation signals and compared to the stand alone GPS approach. without any error control) are depicted in Fig. 5.9(a), whereas Fig. 5.9(b) shows results obtained with the proposed detection/identification/correction filter. In the first case (Fig. 5.9(a)), the final solution is either biased or not bounded during the intervals where errors are present. However, in the second case (Fig. 5.9(b)) the corrected position estimates are in good agreement with the bound thanks to the enhanced scheme.

Experimental data

The previous subsection proved the relevance of the proposed algorithm under totally known conditions (i.e. the exact vehicle position and the magnitude/time of occurrence of interferences were known). This subsection studies the reliability strategy on real experimental data. Its performances in terms of position error and correct bounding will be analyzed. Moreover, dead reckoning data resulting from the WSSs and a gyroscope will be used to show the performance gain obtained when the system is augmented with a second navigation reference. The experimental data was obtained from the same test field campaign described in chapter 4, section 4.8. Indeed, this will enable a performance comparison between the smooth correction criterion proposed in chapter 4 and the de- 5.2 summarizes the performances for all these cases. Analyzing results we can observe that:

• The proposed enhanced EKF/DIC strategy clearly outperforms the standard solution by eliminating aberrant errors and providing an appropriated bounded solution. In this way, the DIC strategy proves to be well-suitable to fight against real urban phenomena.

• As expected, the combination of DR data with the DIC approach provides excellent results. It successfully compensates most of the inconsistent errors in Fig. 5.11(a) and provides bounded solutions. Moreover, when the error correction function is active (see green bars in Fig. 5.11(a)), the position error is reduced. The slight difference between the bounding performance of the GPS and GPS/WSS/Yaw approach might be due to an actual underestimation of the noise values associated to the sensors.

• From the comparison between error tables 4.2-4.3 and 5.2 we can conclude that though both the smooth correction criterion and the DIC strategy are interesting solutions to urban navigation challenges, the last one provides better results. Contrasting performances for the robust GPS/WSS/Yaw hybrid system, the percentage of bounded errors is the same for both strategies but the final mean positioning error is reduced but almost one meter when using the DIC (6.18m vs 6.91m). Nevertheless, it must be noticed that the trade-off for this enhancement in accuracy is a higher computational cost (the smooth correction criterion is easier to implement than the DIC strategy). the presence or absence of line of sight (referred to as LOS and NLOS situations) over the multiple GPS satellites. Therefore two kinds of errors were potentially "corrupting" the pseudoranges, modeled as noise variance or mean value jumps. The time of occurrence and magnitude of these errors were estimated. In this way, realistic measurement models could be obtained. A multiple model EKF was considered as the best adapted solution for this fast-decision/on-line application. Furthermore, the reliability strategy was adapted to exploit the "urban phenomena immunity" of dead reckoning sensors. Simulated and real data validated the relevance of the proposed algorithm showing promising results.

Conclusions and perspectives

The objective of this thesis is to present a reliable 2D GPS-based navigation system particularly adapted to urban environments. The degraded quality of the received signal presents a main challenge when accurate and reliable positioning solutions are demanded. In urban canyon scenarios two main problems are to be considered: a partial or total outage in the satellite visibility, and the lack of integrity in the received pseudo-range measurement. The first problem is addressed by complementing the satellite based navigation system with a dead-reckoning approach. The principle is to use advantageously the redundancy of measurements. The second and most challenging problem is caused by the presence of multipath . Two different approaches are studied with this purpose: firstly by exploiting different types of GPS measurements (such as pseudorange and Doppler measurements) according to their robustness to urban phenomena, secondly by studying and modeling the source of error affecting the commonly used pseudorange signal, so that a consistent noise model can be computed.

Three different hybrid navigation systems have been studied in this thesis. Focusing on land vehicle navigation, GPS augmentations were achieved by the incorporation of inertial sensors (i.e. accelerometers and gyros) and/or odometric data. In particular, on-board wheel speed sensors are of particular interest because they are components of the Antilock Braking System (ABS), and they can be exploited in new generation vehicles at no additional cost. Moreover, as velocity information is obtained for the four vehicle wheels, a self-contained navigation system (based on differential odometry) is possible. A high level integration strategy was proposed to couple GPS, an inertial navigation system (INS) and differential odometric data (referred as WSS). Though the GPS/INS/WSS system dedicated management of the measurements used to compute the solution, that considers the quality of the received GPS signal, was studied in this thesis. Under good reception conditions pseudorange measurements are used. If this is no the case, Doppler measurement are incorporated to the filtering solution in an innovative way. In particular, a weighted Doppler smoothing was proposed so the absolute pseudorange measurements could be merged with the robust Doppler measurements. This approach, referred as smooth correction criterion, presented very good performances when tested with real data. Indeed, it enabled inconsistent navigation solutions to be eliminated and enhancements over the total horizontal position error to be obtained.

The final approach proposed in this thesis doesn't look for an alternative to erroneous pseudorange measurements but it aims at detecting, identificating and correcting the source of error itself. The specific case of multipath interference was considered. Multipath signals can arrive to the receiver either in a line-of-sight (LOS) situation where the direct path is present (i.e. direct visibility over the corresponding satellite), either in a non-line-of-sight (NLOS) situation where the received signal contains reflected components only. We proposed to model the first case by a noise variance jump while the second case was considered to introduce a mean value jump in the noise model. A three hypothesis hierarchical technique was implemented. A first detection of an anomaly presence is followed by the error identification and evaluation. The relevance of the algorithm was validated using real and simulated data. Very promising results in terms of accuracy and correct error bounding were obtained.

Several perspectives can be imagined for this work. In particular the following ones present a special interest:

• Study the benefit brought by additional on-board sensors that can be find in new generation vehicles. In particular, wheel pressure sensors should allow more precise WSS-based velocity measurements while steering sensors should help to determine the heading angle.

• A formalization of the linear merging model between the pseudo-range and Doppler measurements applied in the smooth correction criterion should be investigated. With this purpose, several test field campaign should be done to enable an exhaustive study of the impact urban

List of Publications