Tatiana Sokoltsova

Eric Esbelin
email: eric.esbelin@cea.fr

Marie-Christine Lépy

Quantitative analysis with an energy dispersive X-ray fluorescence instrument equipped with a HOPG filter

Keywords: X-ray fluorescence, quantitative analysis, fundamental parameters, HOPG optic, ray-tracing

We present an energy-dispersive X-ray fluorescence (EDXRF) instrumentation mainly intended for the quantitative analysis of actinides based on their intense L X-ray lines. The experimental setup is equipped with a highly oriented pyrolytic graphite (HOPG) optic as a bandpass filter placed between a sample and a detection system. In such an arrangement, it modifies the spectral distribution of the fluorescence radiation and allows recording of a spectrum in the energy range from 10 keV to 18 keV, while suppressing parts considered to be useless. In order to perform accurate quantitative analysis of the recorded spectra, the classical quantification algorithm based on the fundamental parameters must be adapted to take into account the effect of the HOPG optic on the measured intensities. This requires to determine the transmission function of the HOPG filter and implement it in the algorithm. The validity of this approach is checked in different application examples.

INTRODUCTION

ATALANTE (ATelier Alpha et Laboratoires pour ANalyses, Transuraniens et Études de retraitement) is a nuclear facility of the French alternative energies and Atomic Energy Commission (CEA). Located on the Marcoule nuclear site, it is dedicated to research into spent nuclear fuel reprocessing and management of high-level long-lived radioactive wastes. The ATALANTE analysis laboratory focuses particularly on the analysis of actinides. Measurements are conducted in glove boxes and in hot cells on samples with radionuclides of intermediate and high levels of radioactivity, respectively. An energy dispersive X-ray fluorescence (EDXRF) system was designed and installed in the hot cell of the laboratory for the explicit purpose of performing analyses of samples containing mainly uranium (U) and plutonium (Pu) based on their intense L X-ray lines [START_REF] Esbelin | AIEA Symposium on International Safeguards: Linking Strategy, Implementation and People[END_REF], [START_REF] Esbelin | [END_REF]. The system includes: i) an Rh-anode X-ray tube (4 kW), ii) a cylindrical highly oriented pyrolytic graphite (HOPG) optic, and iii) a high purity germanium (HPGe) detector. HOPG is a mosaic crystal which consists of a large number of small crystallites and owns a high integrated reflectivity due to a large mosaicity [3], [4], [START_REF] Legall | Proceedings of FEL[END_REF]. It is widely used as focusing and dispersion device [START_REF] Kolmogorov | [END_REF], [START_REF]Handbook of Practical X-Ray Fluorescence Analysis[END_REF]. In the present experimental ("nuclearized") setup, the HOPG optic is employed as a bandpass filter between the sample and the detector for specific detection of the fluorescence radiation covering the energy range from 12 keV to 17 keV. The X-rays which are outside of the energy range of interest (e.g., characteristic lines of shielding elements and scattered X-ray tube radiation) are thus suppressed that reduces the emission rate [START_REF] Kolmogorov | [END_REF]. This is especially advantageous for the detection of small contents of analytes which requires to increase the tube current to get results within a short period of time (typically 800 s). In the present setup, the optic element allows selecting the energy region of L Xray lines of actinides such as Th, U, Np, Pu, Am, Cm, but also of K lines of Rb, Pb, Sr, Y, Zr. The system is lead-shielded, since the measurements are carried out on high-activity samples, inducing very high dose rates. The X-ray spectra obtained with the nuclearized setup are quantified with a processing software developed in the laboratory [START_REF] Esbelin | [END_REF] which includes a library of standard spectra. Calibration is performed for each standard element, establishing a nonlinear relationship between the peak intensity and the element concentration [START_REF] Rousseau | [END_REF]. However, the calibration is time-consuming and has to be performed regularly in order to ensure its validity. A spectrum measured with a conventional EDXRF setup can also be processed using an algorithm based on the fundamental parameters (FP) method [START_REF]Handbook of Practical X-Ray Fluorescence Analysis[END_REF], [START_REF] Sharma | X-Ray Spectroscopy[END_REF]. This requires the knowledge of the geometry of the experimental setup, flux and spectral distribution of the excitation radiation, sample matrix, detection efficiency, etc. However, the spectra obtained on the nuclearized setup cannot be quantified in the same manner because the modification of the fluorescence intensities by the HOPG optic is not taken into account in the algorithm. Thus, it is necessary to develop a specific algorithm for the quantitative analysis without standards. This requires to characterize the HOPG response in order to establish a better understanding of its diffraction behavior and define the transmission function. To determine an exact influence of the HOPG filter on the spectral distribution, ray-tracing simulations can be carried out. For the calculations, the exact geometrical parameters of the experimental setup are necessary. Unfortunately, some parameters are not known and the setup cannot be disassembled for such studies for an indefinite period. Thus, this motivated the construction of a miniaturized copy of the nuclearized setup to develop such a FP-based quantification approach and be able to perform tests to get experimental validation. This article aims to discuss the main characteristics of the miniaturized XRF setup equipped with a HOPG optic. The instrument covers the energy range from 10 keV to 18 keV. During this work, the setup was tested in a non-nuclear laboratory using samples with medium-Z elements (Se, Rb, Sr, Y, etc.) whose K X-ray emission lines are in the same energy range as the actinide L X-ray emission lines. A large part of this paper is devoted to the presentation of the optical model which mimics the detection line (from the sample to the detector) of the miniaturized setup using a ray-tracing algorithm. This is applied to calculate the HOPG transmission function which is included in the FP algorithm enabling an accurate quantification. Selected examples are presented to validate this approach.

EXPERIMENTAL SETUP

The studied XRF setup is custom-built and can be considered as a miniaturised copy of the nuclearized system. The principal components of the setup are i) an Ag-anode X-ray tube with a power of 4 W, ii) a HOPG cylinder (Optigraph, Germany) and ii) a silicon drift detector (SDD) (Amptek, USA). The picture of the experimental setup is shown in Figure 1. Picture of the miniaturized XRF setup Transmission-anode X-ray tube with a maximum power of 4 W is used to irradiate the sample. The instrument is positioned at an angle of 22° to the normal of the sample holder at a distance of 16 mm from it. The Bragg diffracted fluorescence radiation is focused by the HOPG optic onto the SDD. The SDD is placed on the axis of HOPG cylinder, which is in the plane perpendicular to that of the X-ray tube and tilted at an angle of 46° with respect to the sample holder axis (see Figure 2). The detection line is based on the von Hamos full-cylinder geometry [START_REF] Anklamm | [END_REF] where the sample is regarded as a geometrical source.

FIGURE 2

Schematic drawing of the experimental setup

To provide radiation shielding, part of the instrument is enclosed in a stainless-steel housing, with a 3 mm thick aluminium (Al) insert to minimize the contribution of fluorescence radiation from the stainless-steel elements. Both the HOPG filter and the silicon detector can be moved separately along their common axis. The samples are placed in a sample cup, which in turn is positioned in an aluminium sample holder. During this work, all samples were prepared from certified standards (1 000 mg L -1 in 2 %HNO3). The HOPG optic is a key element of the experimental setup. It has an inner radius R of 10 mm and a length L is 40 mm. The inner surface of the cylinder is coated with a 200 µm thick HOPG layer with a mosaicity of 0.4° (see Figure 3). The HOPG element was provided with an internal beam stop to avoid the direct transmission of the fluorescence radiation from the sample to the detector without being diffracted [11], [12]. The beam stop is covered with a thin HOPG layer as well. The energy band reflected by the HOPG optic depends on the size of the geometrical source viewed by the HOPG, the length of the cylinder, the sample-to-HOPG and HOPG-to-detector distances and the detection area. With the experimental setup, the fluorescence radiation is collected in the energy range from 10 keV and 18 keV.

HOPG optic

A series of measurement were performed in order to investigate the performances the HOPG filter and demonstrate how this optical element modifies the spectral distribution. Indeed, the HOPG optic can be easily removed to enable recording of the direct X-ray spectrum in the energy range from 3 keV to 30 keV.

FIGURE 4

Fluorescence spectra of a sample containing Se, Rb and Sr at concentrations of 100, 50 and 50 mg L -1 , respectively in HNO3 0.5 M measured without (upper panel) and with (lower panel) HOPG optic

As an example, fluorescence spectra of a sample with Se, Rb and Sr recorded without and with HOPG are compared in Figure 4. With the HOPG, the spectrum in the energy band covered by the filter is intense whereas a low-energy part of the spectrum as well as the coherently and incoherently scattered excitation radiation are considerably reduced.

QUANTIFICATION ALGORITHM

For the quantification of the spectra recorded without HOPG optic, we used the FP-based software PyMCA developed at the European Synchrotron Research Facility (ESRF) [13]. This uses spectrum evaluation presented by Van Espen and Janssens [START_REF] Van Grieken | Handbook of X-ray Spectrometry Second edition[END_REF]. In fact, the fluorescence lines of an element can be grouped (e.g. 𝐾 𝛼 , 𝐾 𝛽 groups, or the entire family of lines (K, L, M etc.) of a given element), and characterised by the total number of counts of all lines, A, in a relevant group. The spectrum of an entire element can be expressed as:

𝑦 𝑝 (𝑖) = 𝐴 ∑ 𝑅 𝑗 𝐺(𝑖, 𝐸 𝑗 ) 𝑁 𝑝 𝑗=1 (1)
where G is the function for a line of energy 𝐸 𝑗 with relative intensity 𝑅 𝑗 . 𝑁 𝑝 represents all the lines in a respective group. The modelling of an entire element permits a reduction in the number of fitting parameters. In the present work, the characteristic X-ray peaks were approximated by a Gaussian function. Other functions such as Lorentzian and Hypermet can also be used to fit the various peaks. The relative transition probabilities of the fluorescence lines between (sub)shells are constant. However, their apparent intensities depend on the absorption in the sample and on absorbers between the sample surface and the active area of the detector. Thus, the attenuation of X-ray must be considered in Eq. ( 1) to define the correct fitting function. In such a case, the apparent intensities are multiplied by the transition probabilities with an absorption correction term. In the case of a polychromatic excitation (e.g., the X-ray tube excitation), the approach considers an incident beam of X-rays of energies 𝐸 0𝑘 with relative rates 𝑤 𝑘 . The relative intensity ratios are defined as:

𝑅 ′ 𝑗 = ∑ 𝑤 𝑘 𝑇 𝐴 (𝐸 0𝑘 )𝜔 𝑗 𝑃 𝑗 (𝐸 0𝑘 )𝑅 𝑗 𝑇 𝐴 (𝐸 𝑗 )[1 -𝑇 𝐷 (𝐸 𝑗 )] 𝜇 𝑇 (𝐸 0𝑘 ) + 𝜇 𝑇 (𝐸 𝑗 ) 𝑠𝑖𝑛𝜑 1 𝑠𝑖𝑛𝜑 2 𝑘 × [1 -𝑒𝑥𝑝 {[ 𝜇 𝑇 (𝐸 0𝑘 ) 𝑠𝑖𝑛𝜑 1 ] + [ 𝜇 𝑇 (𝐸 𝑗 ) 𝑠𝑖𝑛𝜑 2 ] 𝜌𝑑}] (2)
where 𝜔 𝑗 is the fluorescence yield of the given shell, 𝑃 𝑗 considers the possibility of leaving a vacancy in the j shell, 𝜇 𝑇 is the total mass attenuation coefficient of the sample for the primary radiation 𝐸 0𝑘 and the fluorescence radiation 𝐸 𝑗 , 𝜑 1 and 𝜑 2 are the incoming and outgoing angles of the beam with respect to the sample surface, 𝜌 is the sample density, and 𝑑 is the sample thickness. The absorption correction terms 𝑇 𝐴 (𝐸 𝑂𝑘 ) and 𝑇 𝐴 (𝐸 𝑗 ) consider all attenuators between the sample and the active area of the detector for incident and fluorescence energies, respectively. The second absorption correction term [1 -𝑇 𝐷 (𝐸 𝑗 )] represents the absorption of the detector for X-rays. In Eq. ( 2), this apparent transition has to be taken into account with 𝑅 𝑗 ′′ =𝑅 𝑗 ′ / ∑ 𝑅 𝑗 ′ 𝑗 . The total count rate A in an X-ray group of an element is determined by the incident beam intensity 𝐼 0 , the mass fraction C of this element in the sample, and the geometrical efficiency of the detector. The count rate for the group of an element is expressed as:

𝐴 = 𝐼 0 𝐶 𝛺 4𝜋 ∑ 𝑅 𝑗 ′ 𝑗 ( 3
)
where 𝛺 is the solid angle of detection, which is determined by the sample-to-detector distance and the sensitive area of the detector. In Eq.( 3), the total count rate A is a measured value. The flux of the excitation radiation can be determined from the reference material and the geometrical efficiency of the detector can be calculated within the software since all setup parameters are known. In the present experimental setup, in addition, the measured fluorescence intensities are also dependent on the transmission of X-rays through the HOPG optic. Thus, this must be included in Eq. ( 3) as an extra transmission factor.

TRANSMISSION FUNCTION OF THE HOPG FILTER

The investigation of the reflection properties of the HOPG element was performed with the help of the ray-tracing package XRT [START_REF] Klementiev | [END_REF]. The developed optical model includes a source of rays, a HOPG cylinder and an image plane. The sample, being exposed to the primary X-rays, is regarded as the physical source of rays in the simulation model. For the sake of simplicity, the source was modelled as a disc perpendicular to the principal axis of the optical system and parallel to the image plane (see Figure 5). The aperture between the source and the HOPG optical element was created according to the mechanical constructions of the setup which limits the propagation of rays.

Illustration of the ray path from the extended source of radius 𝑟 𝑠 toward the HOPG cylinder passing the circular aperture of radius 𝑟 𝑎 The monoenergetic rays, emitted in a conic ring from a point source lying on the cylinder axis, are reflected from a small crystal surface (blue lines in Figure 5). In the case of an extended source 𝑟 𝑠 , monoenergetic rays can be emitted from any point of the source and consequently with different directions they impinge on the crystal surface at the same axial position. In this case, the angular acceptance of the optical surface is dependent on the size of the source viewed by the crystal, the size of the crystal, and its distance from the source D1. Thus, all effective angles need to be considered in this respect. Rays are incident on the optical surface with an angle θ lying in the range:

𝑎𝑟𝑐𝑡𝑎𝑛 ( 𝑅 -∆𝑟 𝑠 D 1 + 𝐿 ) ≤ 𝜃 ≤ 𝑎𝑟𝑐𝑡𝑎𝑛 ( 𝑅 + ∆𝑟 𝑠 D 1 ) (4)
The effective incident angles dictated by Eq. ( 4) are further limited by the circular aperture of a radius ra at a distance d from the source. Thus, the condition (4) must satisfy the next relation:

𝑎𝑟𝑐𝑡𝑎𝑛 ( 𝑟 𝑎 -∆𝑟 𝑠 𝑑 ) ≤ 𝜃 ≤ 𝑎𝑟𝑐𝑡𝑎𝑛 ( 𝑟 𝑎 + ∆𝑟 𝑠 𝑑 ) (5)
From Eq. ( 4) and Eq. ( 5), the maximum 𝜃 𝑚𝑎𝑥 angle can be calculated to define the angle of the truncated cone delimiting the beam source (see Figure 5, red shaded area). Thus, all possible angles on the crystal surface are taken into account.

Since rays from an extended source propagate with a uniform distribution and with all their directions from the starting point to the optical surface, it is necessary to include a beam stop to prevent direct passing of rays. The exact shape of the beam stop is complicated to reproduce in the simulation model. It can be represented as simple obstacles in the system: two identical circular beam stops on each side of the cylinder and a larger one in the centre, as depicted in Figure 6. The energy window of the collected X-rays is dependent on the solid angle of detection subtended by the detector from the HOPG cylinder, i.e., the detector size and crystal-to-detector distance D2. In order to model the detection area of the SDD, a circular aperture of a size equivalent to that of Si was positioned just in front of the image plane. Thus, the beam collected on the image plane is representative of the detector view.

FIGURE 6

Illustration of the complex optical system for ray -tracing simulation, including all principal components

Let us define the central energy as the energy of the rays reflected from the centre of the HOPG cylinder in the case of a point source (see Figure 5, blue lines). In the geometrical configuration with equal distances D1 and D2, the reflectivity from the first and second order reflections can be calculated. For example, in the case of a point source and a mosaic crystal [START_REF] Sanchez Del Rio | Proc. SPIE[END_REF] to maximize the first order of reflection of 13 keV rays the distances D1 and D2 must be set at 50 mm. In this arrangement, the rays of twice the central energy are reflected in the second order. These rays with higher energies can be reflected from deeper crystal planes, contributing to the output results. Since the experimental measurements were performed at up to over 30 keV, both the first and second orders of reflection had to be computed.

FIGURE 7

Simulation of the reflectivity for the polychromatic beam in the energy range from 5 keV to 35 keV in the first (upper panel) and second (lower panel) orders of reflection. Both simulations were performed for th e central energy of 13 keV in the first order of reflection From Figure 7, it can be seen that rays reflected in the first order of reflection are collected in the energy range from 10 keV to 17 keV, while the rays with energies from 15 keV to 33 keV correspond to the second order of reflection. It is important to point out that almost all optical surface contributes to the first order of reflection, whereas only a narrower part of the crystal reflects the rays in the second order. The full width at half maximum (FWHM) of profiles are 4 keV and 8 keV in the first and second order of reflections, respectively. The intensities of the reflected rays in the second order are much lower, however, and they appear in the energy range of the coherent and incoherent scattering of the excitation radiation. The calculations were performed separately for the same number of rays. The reflection profile of the optical element with the source sampled with a uniform ray distribution was first calculated. At this point, it is necessary to represent mathematically the output of ray distribution after the crystal, in order to integrate it into the PyMCA code. However, since the broadening of the reflection profile varies with the HOPG filter position, the response had to be calculated for each desirable central energy. The reflection profiles of the HOPG optic in the first and second orders of reflection are asymmetric. The complex profile could be fitted with several Gaussian functions [START_REF] Matthew Newville | Release Notes[END_REF], but this would involve a large number of parameters. To simplify the fitting of the filter output and perform it independently of the cylinder position, a spline function was used. This estimates a spline representation of the curve and can be evaluated for a desired set of points. Figure 8 illustrates the complex reflection profile calculated with XRT (blue line) and fitted by spline function (red line).

FIGURE 8

Fit of the reflection profile of the HOPG filter calculated with the source sampled with a uniform ray distribution by means of the ray -tracing method. Calculations were performed in the first and second orders of reflection

RESULTS AND DISCUSSION

Without considering of the transmission function of the HOPG optic, a spectrum is processed incorrectly by the PyMCA package. In Figure 9, it can be seen that the Sr peaks are fitted inaccurately, especially the Sr K-M peak which is largely overestimated by the software. The latter peak is located on the shoulder of the hump of the spectrum, where the transmission is reduced.

FIGURE 9

Processing of the spectrum from the sample with 100 mg L -1 of Sr in HNO3 0.5 M. The spectrum analysis model does not take into account the modification of spectra l distribution by the HOPG filter

In fact, the model can consider the presence of the bandpass filter as a new attenuator between the sample and the detector. To take this into account in Eq. ( 3), each relative intensity 𝑅′ 𝑗 of the line of energy 𝐸 𝑗 is multiplied by a new attenuator term 𝑇 𝑀 𝑗 , which is the relevant transmission function estimated for given positions of the HOPG cylinder and the SDD. Thus, the apparent intensities become:

𝑅 𝑗 * =[𝑅 𝑗 ′ * 𝑇 𝑀 𝑗 ] (6)
This approach was implemented in the PyMca code.

In order to validate the new spectrum analysis model, several experiments were performed with different samples and at various distances D1 and D2. A number of transmission functions for central energies in the range from 12.8 keV to 15 keV were calculated beforehand, to be applied in the model. Figure 10 shows the spectrum acquired for the solution containing 100 mg L -1 of Se in HNO3 0.5 M (black line) and the result of the fit (red line). The comparisons of the prepared and estimated mass fractions are reported in Table 1 (Processing №1). To go further in the analyses, the spectrum from a sample with three elements at different concentrations was processed (see Figure 12). It can be seen that the lines on the shoulders of the hump, where the reflectivity is decreased, fit perfectly. The concentrations of all elements in the sample could be accurately calculated, despite intensities reduced by the filter (see Table 1, Processing №3).

FIGURE 12 Processing of the spectrum from a sample with 100, 50, and 50 mg L -1 of Se, Rb, and Sr, respectively, in HNO3 0.5 M

CONCLUSIONS

We demonstrated the possibility of analysing spectra obtained with an EDXRF setup equipped with HOPG optic, using a quantification algorithm based on the fundamental parameters. The HOPG optic device was used as a broadband filter, which favoured the collection of the fluorescence radiation in the desired energy ranges, and suppressed the low energy peaks and scattered radiation of the excitation source. The selected energy domain extended from 10 keV to 18 keV, which corresponds to that of the actinide L X-ray lines.

For the accurate quantitative analysis of the recorded spectra, the transmission function of the HOPG cylinder (taken into account both the first and second orders of reflexion) was calculated by means of the ray-tracing package XRT. This enabled us to represent the entire detection channel of the experimental setup. The calculated transmission functions, fitted with spline functions, were implemented in the FP-based code PyMCA in such a way that each intensity of each measured line of a given energy was multiplied by the relevant transition efficiency. The new model was tested for different concentrations of medium-Z elements (Se, Rb, Sr, Y). While the direct quantification, i.e. without the calculation of HOPG transmission, was strongly biased with relative differences up to 50 % compared to the prepared mass fractions. Whereas, the new approach successfully estimated mass fractions with a relative difference less than 5 % from the prepared values. More details on this work can be found in the exhaustive study [START_REF] Sokoltsova | [END_REF]. In the future, this approach could be applied to the

FIGURE 1

 1 FIGURE 1

FIGURE 3

 3 FIGURE 3

FIGURE 10

 10 FIGURE[START_REF] Anklamm | [END_REF] Processing of the spectrum from the sample containing Se at a concentration of 100 mg L -1 in HNO3 0.5 M

FIGURE 11

 11 FIGURE11 Processing of the spectrum from the sample containing Rb and Y at concentrations of 100 mg L -1 in HNO3 0.5 M

Table

ACKNOWLEDGEMENTS

We are grateful to Konstantin Klementiev from MAX IV Laboratory (Lund, Sweden) for his help and valuable advices. We acknowledge also the team of Prof. Birgit Kanngießer from the Berlin Laboratory for innovative X-ray technologies (BLIX) and notably Richard Gnewkow for fruitful discussions on the application of HOPG optics.

nuclearized instrumentation in order to avoid time-consuming calibration steps.


Joelle Al Hage

Philippe Xu

Philippe Bonnifait

Javier Ibanez-Guzman

Localization Integrity for Intelligent Vehicles through Fault Detection and Position Error Characterization

Keywords: fault detection and exclusion, error characterization, protection level computation, map aided localization

Localization integrity consists in providing a realtime measure of the level of trust to be placed in the localization estimates as vehicles operate. It provides a means of knowing whether position estimates are usable for navigation purposes. This paper formalizes the integrity concept and its underlying principles. Vehicles operate in different navigation environments, and so multiple sensors are used to ensure the required performance. Different sources of error exist. They must be bounded according to the acceptable level of risk for the application.

This paper presents a generic approach for addressing integrity. It combines measurement rejection (for measurements considered to be faults) and position error characterization. For this purpose, a multi-sensor data fusion with a Fault Detection and Exclusion algorithm is constituted using a bank of information filters. These filters allow detected faults to be isolated without any prior assumption regarding the number of simultaneous errors. In addition, external integrity is expressed as a Protection Level of the localization solution. It uses a Student's t-distribution in order to bound the distribution of the position error applicable to small integrity risks after a learning step. The approach is tested on data acquired on public roads using an experimental vehicle equipped with off-the-shelf proprioceptive and exteroceptive sensors together with an HD map. The results obtained validate the proposed approach.

I. INTRODUCTION

For mobile platforms, knowing their location is essential for navigation purposes. To increase the accuracy of localization solutions, multi-sensor data fusion methods provide good results. Higher levels of vehicle autonomy not only need higher accuracy, but also an ability to estimate error bounds in order to address uncertainty in the localization estimates in real time. The integrity concept includes this second feature. It is defined as the measure of trust on the correctness of the location estimates applied to navigation systems [START_REF] Ober | Integrity prediction and monitoring of navigation systems[END_REF]- [START_REF] Zabalegui | A review of the evolution of the integrity methods applied in gnss[END_REF].

The accuracy of GNSS-based absolute positioning systems is improved through the fusion with multiple sources of information. For road vehicles, their location is established with respect to a map that incorporates a stored geometric description of the roads [START_REF] Betaille | Creating enhanced maps for lanelevel vehicle navigation[END_REF]. The estimated positions are associated with road maps via the map matching process [START_REF] Quddus | Current mapmatching algorithms for transport applications: State-of-the art and future research directions[END_REF], [START_REF] Kasmi | Map matching and lanes number estimation with openstreetmap[END_REF]. Information from exteroceptive sensors (e.g., cameras) is used to detect navigation feature landmarks. These are georeferenced and stored in the map to enhance the positioning output [START_REF] Tao | Mapping and localization using GPS, lane markings and proprioceptive sensors[END_REF]. 1 Heudiasyc, UMR 7253, Université de Technologie de Compiègne, France. 2 Renault S.A. 1 Av. Du Golf, 78064 Guyancourt, France.

Data fusion is usually done using a Kalman Filter (KF). The estimated states are computed along with a covariance matrix that encodes the inherent uncertainty of the estimates. This matrix statistically evaluates the estimation error.

The accepted risk level can be defined by the Target Integrity Risk (T IR), which depends on the navigation context and any applied safety policy. The Protection Level (PL) is defined as the statistical bound on the Position Error (PE) such that the probability of the absolute error exceeding the PL (empirical integrity risk denoted IR) is smaller or equal to the T IR. Therefore, for a chosen T IR, the system meets the integrity requirement if the IR is lower than the T IR:

Pr (PE > PL) ≤ T IR. (1)
The system is unavailable if PL>AL, where AL is the Alarm Limit defined as the largest error allowable for safe operation without issuing an alert. The different integrity concepts are often shown using a Stanford-ESA diagram. This plots PE against PL for a set of measurements with a given T IR, which is helpful in assessing the performance and usability of a localization system [START_REF] Tossaint | The stanford ESA integrity diagram: A new tool for the user domain SBAS integrity assessment[END_REF].

For autonomous vehicles, it is a question of choosing the correct T IR level. In the ESCAPE project 1 , the integrity levels for Galileo-based GNSS hybrid localization systems were estimated. The T IR chosen was 10 -6 per estimated location. This is for an accident probability of about 10 -9 . Similar figures were chosen in the SaPPART COST Action 2 . The hypothesis applied was that no accident should occur due to localization errors.

In practice, it is often observed that the uncertainty estimates computed from the Kalman filter are suboptimal and not consistent. That is, PLs computed from the covariance matrix do not satisfy equation [START_REF] Ober | Integrity prediction and monitoring of navigation systems[END_REF]. This is because one or more fundamental hypotheses of the Kalman filter framework are violated. The sources of errors include linear approximations, non-Gaussian noise, non-centered noise, correlated errors, outliers, etc.

The present paper proposes a novel localization solution from an integrity perspective, and the bounding of the estimation error. Our proposed approach combines measurement rejection and position error characterization. A Fault Detection and Exclusion (FDE) step is first applied, then a PL is computed using a fine-tuned Student's t-distribution optimized from recorded data.

The main contributions of the paper are i) using a bank of information filters for the data fusion, with an FDE that is able to handle different sensors and residuals in the state space, ii) including a Student's t-distribution for PL computation, chosen according to a learning strategy, iii) combining measurement rejection and position error characterization in an integrity framework without the assumption of one fault at a time.

The paper is organized as follows. Section II presents a state-of-the-art of integrity applied to localization systems. The rationale behind the proposed approach is presented in section III. The multi-sensor data fusion strategy is described in section IV. The integrity study, detailed in Section V, includes the FDE step and the computation of the PL. The approach is tested using data acquired in public roads. The results are included and analyzed in Section VI.

II. RELATED WORK

Integrity determines the level of trust that can be placed in localization estimates. This can be expressed in terms of PLs defined by the acceptable integrity risk, as in equation [START_REF] Ober | Integrity prediction and monitoring of navigation systems[END_REF]. Most studies on integrity focus on the aerospace field. This section looks into the different methods that have been applied to determine levels of integrity. Three types of methods are considered: the rejection of measurements (that is, FDE), error characterization according to the source of error, and mapaided integrity methods applied to road vehicles [START_REF] Cosmen-Schortmann | Integrity in urban and road environments and its use in liability critical applications[END_REF], [START_REF] Worner | Integrity for autonomous driving: A survey[END_REF].

The main drawback of the FDE integrity approach when applied to intelligent vehicles is that more than one measurement may be faulty at the same time, and may therefore fail to provide a navigation solution with an associated PL. Error characterization approaches consist in modeling the measurement error distribution and computing a PL without an FDE step. They can be applied either in the range measurement domain [START_REF] Ober | Overbounding the SBAS integrity equation[END_REF] or in the position domain [START_REF] Ahmad | Characterization of GNSS receiver position errors for user integrity monitoring in urban environments[END_REF]. Error characterization approaches can lead to large PLs in the presence of faulty measurements, increasing the unavailability of the localization system.

A. Classical integrity methods and their limitations for intelligent vehicles

In aviation, PLs are measurement-independent. They depend on satellite geometry and the prior distribution of the errors. These PLs are therefore predictable at any point on Earth [START_REF]Minimum operational performance standards for global positioning system/wide area augmentation system airborne equipment[END_REF]. A classical application of the FDE approach to the GNSS localization problem is Receiver Autonomous Integrity Monitoring (RAIM). This monitors the integrity of the position using redundant information (e.g., pseudo-distances) to detect satellite faults. Conventional RAIM relies on a one-fault-at-atime detection scheme. New RAIM methods incorporate FDE mechanisms that allow for service continuity [START_REF] Nikiforov | Advanced raim algorithms: First results[END_REF]. Integrity monitoring within the RAIM can be achieved through four different methods: least squares residuals [START_REF] Parkinson | Autonomous GPS integrity monitoring using the pseudorange residual[END_REF], parity space [START_REF] Sturza | Navigation system integrity monitoring using redundant measurements[END_REF], range comparison [START_REF] Lee | Analysis of range and position comparison methods as a means to provide GPS integrity in the user receiver[END_REF], and maximum solution separation [START_REF] Brown | Self-contained GPS integrity check using maximum solution separation[END_REF].

The best-known method for PL computation is the FDE slope-based method [START_REF] Brown | A baseline GPS RAIM scheme and a note on the equivalence of three RAIM methods[END_REF]. The test statistic, constructed in the observation space, follows a central Chi-squared distribution in the non-faulty case and a non-central Chi-squared in the faulty case. It assumes a single biased pseudo-range at a time. The threshold is chosen such that a given false alarm probability Pr f a is guaranteed. The PL is obtained by computing the effect of the minimum detectable bias on the horizontal PE.

Another method is Maximum Solution Separation (MSS). This requires a rejection step [START_REF] Brown | Self-contained GPS integrity check using maximum solution separation[END_REF], [START_REF] Brenner | Integrated GPS/inertial fault detection availability[END_REF]. MSS is based on the construction of a test statistic in the position domain. It also assumes one fault at a time.

Expecting one fault at a time is reasonable for aviation applications, but not for ground vehicles [START_REF] Marchand | Characterization of gps multipath for passenger vehicles across urban environments[END_REF]. Where there are two faulty satellites simultaneously, computing the PL requires examining every pair of satellites, which for N satellites gives N! (N-2)!×2! combinations to be tested [START_REF] Liu | Theoretical analysis of RAIM in the occurrence of simultaneous two-satellite faults[END_REF]. When applied to GNSS constellations (e.g. GPS, Galileo, Glonass, etc.), this can become computation-intensive and thus difficult to use.

In [START_REF] Roturier | The SBAS integrity concept standardised by ICAO-application to EGNOS[END_REF], the PL computation is done using the error characterization of two parameters: the User Equivalent Range Error (UERE) and the Dilution Of Precision (DOP). The UERE calculates the quality of the range measurements based on ephemeris uncertainties, receiver noise, propagation errors, and clock and timing errors. It is modeled as a Gaussian random variable with a zero mean and a standard deviation σ UERE . The DOP represents errors that are due to the satellite spatial geometric distribution. The assumption of a zero mean Gaussian distribution for GNSS errors, even when operating with differential corrections in urban environments, corresponds only rarely to reality, due to the presence of multipath and biases. The same applies to Satellite-Based Augmentation Systems (SBAS) like EGNOS [START_REF] Ali | Characterizing local effects on protection level concept in urban environments[END_REF] [25] [START_REF] Marais | Land mobile GNSS availability and multipath evaluation tool[END_REF].

The PLs obtained in aviation applications are too conservative for vehicular applications [START_REF] Cosmen-Schortmann | Integrity in urban and road environments and its use in liability critical applications[END_REF]. Airplanes travel in 14-kilometer-wide air corridors [START_REF]Instrument Procedures Handbook[END_REF]. Simulations shows that while RAIM availability is around 100% in open sky, it decreases dramatically in urban environments because of the occlusion of satellite signals by high-rise buildings, even with a 100-meter-wide alert limit [START_REF] Nowak | Snapshot RAIM algorithms availability in urban areas[END_REF]. RAIM availability is 55% in urban environments when using the GPS-Galileo dual constellation and with Pr f a = 5 × 10 -3 and Pr md = 5 × 10 -5 . With the GPS constellation alone, this availability decreases to 7% [START_REF] Salos | Receiver autonomous integrity monitoring of GNSS signals for electronic toll collection[END_REF].

Localization systems used for intelligent vehicles are enhanced by incorporating other proprioceptive and exteroceptive sensors like video cameras or LiDARs. In urban conditions, measurements are affected by unpredictable errors, and so the computation of the PL must include the residuals on the actual measurements.

Given these limitations, new methods for integrity monitoring in road conditions are needed. The objective is to address the presence of multiple faults simultaneously, the availability of exteroceptive sensors other than GNSS, and the use of maps that store landmarks and features that can be used for localization. Attention should also be paid to non-Gaussian distributions of the errors.

B. Integrity applied to road vehicles

Different integrity estimation methods applicable to urban environments have been proposed. In [START_REF] Margaria | The local integrity approach for urban contexts: Definition and vehicular experimental assessment[END_REF]- [START_REF]Proof-of-concept of the local integrity approach: Prototype implementation and performance assessment in an urban context[END_REF], the idea is to characterize and predict the degradation of the GNSS signals by implementing a collaborative spatio/temporal approach between multiple vehicles within a Vehicular Ad-hoc Network (VANET). The method estimates an effective σ UERE using the covariance of the residual vector on a specific position. It does not assume a generic model, and includes the effects of noise and multipath on the quality of pseudo-range measurements for a specific position. However, it requires multiple vehicles, or multiple runs in the same place.

Another approach characterizes errors in the position domain rather than modeling individual pseudo-range measurement errors [START_REF] Ahmad | Characterization of GNSS receiver position errors for user integrity monitoring in urban environments[END_REF]. The Horizontal Position Error (HPE) is computed from experimental data with respect to a ground truth. The HPE is fitted by a Rayleigh distribution in open sky, using an empirical Probability Density Function (PDF) and a Cumulative Density Function (CDF). In urban environments, the pseudo-range errors are likely to be non-Gaussian due to multipath. This implies that errors generate distributions with heavy tails. To provide flexibility, the authors chose the generalized Pareto distribution. Once the CDF of the HPE is determined, the PL can be obtained from the inverse of the CDF for a chosen missed detection probability.

The characterization of the position errors is also proposed in [START_REF] Madrid | Device and method for computing an error bound of a kalman filter based gnss position solution[END_REF], [START_REF] Pedro | New approach for integrity bounds computation applied to advanced precise positioning applications[END_REF]. The method is known as the Kalman Integrated Protection Level (KIPL). It is an extension of an earlier method, namely the Isotropy-Based Protection Level (IBPL) [START_REF] Saenz | Method for autonomous determination of protection levels for gnss positioning based on navigation residuals and an isotropic confidence ratio[END_REF], [START_REF] Azaola | Autonomous isotropybased integrity using GPS and GLONASS[END_REF]. The KIPL uses a KF for state estimation and error propagation. To calculate the PL, errors are assumed to be distributed according to a multivariate Student's tdistribution. The total error of the KF is considered to be the sum of contributions from several error sources. For each measurement i, the position error is separated into two distinct elements: an error of the measurement i projected to the state model, and an error of the measurement i accumulated over previous epochs. Given that the sum of two t-distributions is not t-distributed, the authors use a fitting method to adjust the sum to a t-distribution. Once the distribution is known, the associated protection level PL i can be computed and the PL is considered to be the sum of PL i . To derive the KIPL, the authors make certain assumptions, such as uncorrelated Kalman gain and a diagonal innovation covariance matrix.

The various methods mentioned above characterize errors without removing large biases, that is to say no FDE is applied. These biases may increase the unavailability of the system.

For the sake of completeness, we also need to mention emergent map-aided integrity methods. Empirical methods that use map-matching to estimate integrity are proposed in [START_REF] Quddus | Integrity of mapmatching algorithms[END_REF]. The integrity of the map matching is considered as high if the uncertainty associated with the position solution is small, and if the absolute differences between the vehicle pose and the road link (in the map) in terms of distance and heading are small. For this purpose, a fuzzy logic model is used to derive a metric between 0 and 100 that reflects the map matching confidence.

[38] used precise point-positioning GNSS and a monocular camera together with an HD map to study the precision and reliability of the vehicle state estimation in the presence of errors. This approach examined internal reliability (FDE) and external reliability (the impact of undetected faults).

Another approach addresses lane level integrity [START_REF] Toledo-Moreo | Lane-level integrity provision for navigation and map matching with GNSS, dead reckoning, and enhanced maps[END_REF]. GNSS and vehicle odometry are merged with a digital map through a Particle Filter (PF). The integrity of the map-matched position is represented by two variables. The Lane Occupancy Probability (µLO) of a segment is obtained by adding the weights of the particles corresponding to that segment. The confidence level of positioning and map matching is estimated. This approach also estimates the Lane Positioning Protection Level (LPPL) defined by kσ pos , where k is calculated using a Rayleigh distribution according to a missed detection probability and σ pos obtained from the maximum eigenvalue of the covariance matrix along the x and y-coordinates. The position is valid if µLO > µLO threshold and LPPL < LPPL threshold .

Integrity monitoring can also be derived using a PF that estimates a set of likely map-matched hypotheses, according to the PL provided by a GNSS receiver [START_REF] Li | Map-Aided Dead-Reckoning With Lane-Level Maps and Integrity Monitoring[END_REF]. That is, the current fix of the GNSS receiver is used to select the coherent hypotheses using a Mahalanobis distance. Where the remaining set contains a single hypothesis (as opposed to multiple hypotheses or no hypotheses at all), it may be decided to use the position estimate.

III. PROPOSED APPROACH

It is our assertion in this paper that localization systems for intelligent vehicles need to use redundant and complementary technologies. GNSS systems are indispensable, given their Earth-wide coverage, their constantly improving positioning quality and the low cost of receivers. These systems are hybridized using inertial systems or Dead-Reckoning (DR) using wheel speed measurements, for instance. HD maps are a valuable complementary source of positioning information. By combining the map with measurements from exteroceptive sensors such as LiDARs and cameras, the car can derive very accurate localization information. This paper formulates an integrity concept applicable to any multi-sensor localization system.

The approach addresses the constraints encountered in land navigation, such as errors due to multipath and NLOS of GNSS satellites [START_REF] Braasch | Multipath effects[END_REF], [START_REF] Groves | Intelligent urban positioning using multi-constellation GNSS with 3d mapping and NLOS signal detection[END_REF], measurement errors of the exteroceptive sensors, landmarks and features errors in the map, and data association issues.

Our proposed approach combines measurement rejection and position error characterization. It merges the information provided by all information sources. The data fusion of a DR system with a low-cost single-frequency GNSS receiver is considered. A smart camera is used to provide measurements of detected lane markings georeferenced in an HD map (see Figure 1).

Measurements are affected by errors mainly due to the navigation environment, although errors on the camera measurements may originate in the measurement itself or in the map. These errors have a direct effect on the integrity of the localization if they are not detected and excluded from the fusion procedure. Therefore, an FDE step is added in order to test the consistency of the exteroceptive measurements with the DR measurements. Our approach is based on the use of the Information Filter (IF) for state estimation, and Mahalanobis distance in the state space for fault detection. The IF is mainly advantageous for multi-sensor data fusion and for FDE as a result of the summation part that appears in the update step [START_REF] Hage | Improved outdoor localization based on weighted kullback-leibler divergence for measurements diagnosis[END_REF].

Once all the detected faulty measurements are removed, an error model is built to compute a consistent PL. For this step, we use a Student's t-distribution. This is characterized by a degree of freedom (dof) that shapes the tail of the distribution, leading to a more consistent PL for small integrity risk (e.g., 10 -3 ). The heavy tail of the Student's t-distribution handles large errors with low occurrence probabilities, and this cannot be done when using a Gaussian distribution. The tuning of the dof can be done experimentally in a learning step according to the navigation environment (urban, suburban, open sky, etc.). This is in line with the strategy for deploying initially autonomous vehicles within a well-defined operational design domain at first. For intelligent vehicles, PLs are expressed in the Along Track (AT) and Cross Track (CT) directions that describe the relative position of the vehicle with respect to the road. This allows different alarm thresholds depending on the direction.

IV. MULTI-SENSOR DATA FUSION

At instant k, the state vector is considered to be the position and the heading of the vehicle in an East-North-Up (ENU) frame denoted R O :

X = [ x y θ ] T (2)
The propagation equation is obtained using DR as follows:

X k+1/k = f (X k/k , u k ) + v k =   x k/k + cos(θ k/k + ω k 2 )∆ k y k/k + sin(θ k/k + ω k 2 )∆ k θ k/k + ω k   + v k (3)
where u k = ∆ k ω k T is the input vector composed of an elementary rotation ω k obtained from the gyro and an elementary displacement ∆ k obtained from the wheel speed sensors, and v k is modeled as Gaussian white noise with covariance matrix Q k .

The information matrix, defined as the inverse of the covariance matrix, is written as:

Y k+1/k = P -1 k+1/k = (F k P k/k F T k + B k (Q u ) k B T k + Q k ) -1 (4)
where Q u is the covariance matrix corresponding to the input vector; F k and B k are the Jacobian matrices computed as

F k = ∂ f ∂ x | X=X k/k and B k = ∂ f ∂ u | u=u k .
The information vector is obtained from the state vector and takes the form:

y k+1/k = Y k+1/k X k+1/k (5)
The information matrix and information vector are updated according to equations ( 6) and [START_REF] Tao | Mapping and localization using GPS, lane markings and proprioceptive sensors[END_REF]:

Y k/k = Y k/k-1 + n ∑ i=1 I i,k (6)
y k/k = y k/k-1 + n ∑ i=1 i i,k (7)
where n is the number of observations and I i,k and i i,k are the information contributions associated with the observation Z i . The computation of these values is detailed in the following. Please note that the IF is mathematically equivalent to the KF. However, the main advantage of IF over KF can be seen in equations ( 6) and [START_REF] Tao | Mapping and localization using GPS, lane markings and proprioceptive sensors[END_REF], where the update step is modeled as a simple summation. This facilitates the FDE step.

The observations from the GNSS receiver and the smart camera are merged with the DR measurements to improve the localization so as to reach a lane-level accuracy.

The fusion of DR with the GNSS is done through a loosely coupled architecture. Therefore, the observation vector corresponding to the GNSS measurement is written as:

Z GNSS = x GNSS y GNSS = t x cos θ -t y sin θ + x t x sin θ + t y cosθ + y , (8)
where t x and t y are the translation of the antenna (frame R G ) with respect to the body frame (R B ) located at the middle of the rear wheel axis of the vehicle where the pose (x, y, θ ) is defined (Figure 2). The GNSS information contribution is given by:

I GNSS,k = H T GNSS,k R -1 GNSS,k H GNSS,k (9)
i GNSS,k = H T GNSS,k R -1 GNSS,k Z GNSS,k -ẐGNSS + H GNSS,k X k/k-1
where H GNSS is the observation matrix:

H GNSS = 1 0 -t x sin θ k/k-1 -t y cos θ k/k-1 0 1 t x cos θ k/k-1 -t y sin θ k/k-1 (10)
and R GNSS is the covariance matrix of the observation noise.

The smart camera detects lane markings and returns the coefficients of a Taylor expansion of a clothoid in a virtual frame R M located at the center of the front bumper, while the camera itself is located behind the windscreen:

y = C 3 x 3 +C 2 x 2 +C 1 x +C 0 ( 11
)
The camera is able to detect up to four markings at a time, two on each side: right and left. In this paper, we consider only the use of the parameter C 0 (the lateral distance between the lane marking and the frame R M ):

Z i = C 0 ( 12
)
where i can represent C 0,left,1 , C 0,left,2 , C 0,right,1 , C 0,right,2 if they exist.

The camera observation model developed in [START_REF] Tao | Mapping and localization using GPS, lane markings and proprioceptive sensors[END_REF] can be written as:

C 0 = (P x sin θ + y -y A ) • x AB -(P x cos θ + x -x A ) • y AB x AB • cos θ + y AB • sin θ (13)
where the lane marking corresponding to the map-matched segment

[AB] at coordinates [x L , y L ] is expressed in frame R O as: O L = O T B • B T M • M L ( 14
)
with O T B and B T M the transformation matrices from R B to R O and from R M to R B , respectively. Therefore,

x L y L = P x cos θ +C 0 sin θ + x P x sin θ -C 0 cos θ + y ( 15
)
where M L = [0,C 0 ] and P x is the distance between R B and R M (Figure 2).

V. INTEGRITY PROCESSES

A. Fault detection and exclusion

The observations from the GNSS receiver or from the camera can be faulty (or erroneous). The origin of the fault can come from the sensor, the map or the environmental conditions. Erroneous measurements (faults), if not detected, directly affect the localization quality. Hence, an FDE step is added, where the Mahalanobis distance between the predicted state and the updated state is used as a residual:

r k = (X k/k -X k/k-1 ) T Y k/k (X k/k -X k/k-1 ) (16)
X k/k and Y k/k are obtained from a principal filter that uses all available measurements at instant k (cf. Figure 3). Note that this residual acts in the state space and detects faults which have a direct influence on the position estimate. Faulty measurements that have a negligible effect on position are not excluded from the fusion procedure. In some cases, using measurements with errors can be preferable to isolating them [START_REF] Betaille | 3d-city-model-aided gnss accurate positioning with integrity provision using simplified geometry of buildings[END_REF].

If r k exceeds a threshold, a fault is detected and its origin has to be determined. The threshold is obtained from a Chisquared distribution with 3-degrees of freedom according to a given false alarm probability (Pr f a =0.05 in this paper).

For the exclusion step, a set of residuals r i is generated using a bank of information filters EIF i :

r i,k = (X k/k,i -X k/k-1 ) T Y k/k,i (X k/k,i -X k/k-1 ) (17)
where X k/k,i and Y k/k,i are obtained from a filter (EIF i ) that uses only the observation i in the update step (i ∈{GNSS, C 0,left,1 , C 0,left,2 , C 0,right,1 , C 0,right,2 }). If r i,k exceeds the threshold value, the measurement i is excluded from the fusion procedure by simply subtracting its information contribution from the principal filter (Eq. ( 6) and ( 7)). For r i,k the threshold is also obtained from a Chi-squared distribution with 3-degrees of freedom. In this paper, for simplification, the false alarm probability is chosen to be the same for the different residuals.

In practice, a threshold optimization has to be done to improve the detection quality.

At a given time the complete set of residuals is not always available, since the sensor measurements are managed asynchronously, but this is not an issue for the method. For instance, the availability of the camera measurements depends on the navigation environment.

Using this approach, a distinction can be made between camera faults (including map-matched errors) and map faults. Map faults refer to errors in the HD map. As the camera gives an indication of the quality of the lane marking detection, a lane marking with a quality below a given value is

DR Sensors

𝑋 𝑘/𝑘 (5) , 𝑌 𝑘/𝑘 (5) rejected. This initial filtering considerably reduces the number of outliers. Given that the camera is able to detect up to four lane markings, two on each side, the redundancy can be used to differentiate between camera and map faults. Let us take the following example for illustration. At instant k, suppose the camera observes simultaneously two lane markings on the left (and so provides two measurements C 0,left,1 and C 0,left,2 ). If the residual r C 0 ,left, 1 exceeds the threshold value, we can conclude that something is wrong with this measurement. At this stage, we cannot specify if the error comes from the camera or from the map. Suppose now that the residual r C 0 ,left, 2 is below the threshold. In this case, at the same instant and on the same side, the camera has one good match and one bad. In these conditions, it is likely that the error comes from the map, since the two measurements are provided with high confidence. This is a phenomenon that we have observed several times by comparing the results using the HD map. A flag can be stored in the map to indicate that the map-matched segment is potentially erroneous. Now, if r C 0 ,left,2 is greater than the threshold, it cannot be concluded either that there is or that there is not an error in the map. Indeed, if the two observations on the same side are faulty, the error may be due to the environment, the camera itself, or the map-matching procedure. Likewise, if only one measurement is available on a given side, it is not possible to differentiate between camera faults and map faults. Finally, if all the residuals r i indicate the presence of errors, the hypothesis of a filtering problem or DR problem is likely and an alarm has to be raised. The origin of the error (camera or map) does not affect the exclusion procedure. If r C 0 , j indicates the presence of a fault, the measurement C 0, j is not used in the update step, irrespective of the origin of the measurement error.

GNSS Receiver

B. Protection level computation

After having removed the detected faulty measurements, the objective is now to bound the estimation errors. The aim is to avoid an underestimation of the PL that could be due to undetected biases, non-Gaussian error distributions, or linearization of non-linear models that could result in an underestimation of the covariance matrix. Hence, when computing the PL, the multivariate Gaussian distribution is replaced by a multivariate Student's t-distribution, which has a heavier tail than the Gaussian distribution and is able to handle errors with small probabilities that are not taken into account when using the Gaussian distribution [START_REF] Dhital | A new approach for improving reliability of personal navigation devices under harsh gnss signal conditions[END_REF].

A random variable X is said to follow a multivariate tdistribution with mean vector µ, shape matrix R S and dof ν, if it can be written in the form [START_REF] Kotz | Multivariate T-Distributions and Their Applications[END_REF]:

X = µ + ω -1 2 y ( 18
)
where ω is distributed according to a Gamma distribution Γ ν 2 , ν 2 and y has a zero-mean Gaussian distribution with covariance P. It is denoted as X ∼ St(µ, R S , ν).

The dof of the Student's t-distribution has a direct effect on the tails of the distribution and if ν → ∞ the t-distribution converges to a Gaussian one. The covariance of the t-distribution is defined only if ν > 2, and in this case it can be written as P = ν ν-2 R S . Figure 4 shows, for a given target integrity risk (T IR = 10 -3 ), the interest of a t-distribution for PL computation as compared to a Gaussian distribution. The blue plot represents a t-distribution with a shape matrix equal to the covariance of the Gaussian distribution. Here the PL obtained from this t-distribution is very conservative in comparison to that obtained from the Gaussian. The green plot represents a tdistribution with covariance equal to the covariance of the Gaussian distribution. For small integrity risk, the PL obtained from this t-distribution is more conservative than that obtained using the Gaussian distribution, but it is less conservative than that illustrated in blue. For intelligent vehicles, the aim is to avoid an undervaluation of the PL under small integrity risks, without however being too conservative; otherwise the system will become unavailable. For this purpose, the matrix P k/k will be considered as the covariance matrix of a multivariate Student's t-distribution with a fixed dof ν.

The horizontal protection level is expressed as follows [START_REF] Madrid | Device and method for computing an error bound of a kalman filter based gnss position solution[END_REF]:

PL k (α) = K(α, ν) ν × max(eigenvalue(R S )) (19)
where R S is the shape matrix of the t-distribution with a dof ν. Given that the covariance can be written as P k/k = ν ν-2 R S , the PL can be expressed as:

PL k (α) = K(α, ν) √ ν -2 max(eigenvalue(P k/k )), ( 20
)
where K is obtained from a multivariate t-distribution with dof ν and according to a given confidence level α. The K(α, ν) value can be computed by numerically solving the following equation [START_REF] Madrid | Device and method for computing an error bound of a kalman filter based gnss position solution[END_REF]:

α = 2 B d 2 , ν 2 ∞ K y d-1 (1 + y 2 ) ν+d 2 dy. ( 21
)
where B d 2 , ν 2 is the Beta function and d is the dimension of the t-distribution; for horizontal PL, d = 2.

For intelligent vehicles, we are interested in PLs in the AT and CT directions. Therefore, the covariance matrix in the ENU frame is projected to the (AT, CT) frame. [START_REF] Ober | Integrity prediction and monitoring of navigation systems[END_REF][START_REF] Zhu | GNSS position integrity in urban environments: A review of literature[END_REF] ) be the eigen values and eigen vectors of the 2 × 2 position component of the covariance matrix P k/k . The expression of V 1 and V 2 in the (AT, CT) frame is given as:

Let (λ 1 , λ 2 ) = eigenvalue(P k/k(1,2),(1,2) ) and (V 1 ,V 2 ) = eigenvector(P k/k(1,2),
V 1,AT V 1,CT = R pro j •V 1 and V 2,AT V 2,CT =R pro j •V 2
where

R pro j = cos θ sin θ -sin θ cos θ (22)
with θ the orientation of the vehicle.

The AT and CT PLs are obtained as:

PL AT = K(α, ν) √ ν -2 max (|λ 1 V 1,AT | , |λ 2 V 2,AT |) ( 23
)
PL CT = K(α, ν) √ ν -2 max (|λ 1 V 1,CT | , |λ 2 V 2,CT |) (24)
The errors in the AT and CT directions are expressed in terms of the error in the East (er x ) and North (er y ) directions as:

er AT er CT = R pro j er x er y

The dof of the Student's t-distribution is defined based on a learning strategy which includes prior knowledge of the driving context. Further details are provided in the next section.

VI. EXPERIMENTAL STUDY

Experiments were performed in the city of Rambouillet (figure 5). The test vehicle is equipped with a Ublox 8T single frequency multi-constellation GNSS receiver and a production camera for lane detection (Mobileye EQ3). An HD-map of the area, designed for testing autonomous vehicles, is also available. The location ground truth is generated using a high-end system (NovAtel SPAN-CPT) with RTK corrections. Acquisitions were made at 50 Hz for the ground truth, 50 Hz for the dead-reckoning, 2 Hz for the GNSS fix, and 35 Hz for the camera (sub-sampled at 3.5 Hz to reduce time correlation errors).

Three trajectories with several runs were carried out and separated into training and testing trajectories. Each of the three trajectories was repeated twice, the first run for the learning phase (training trajectory) and the second for the test phase (test trajectory). Trajectory #1 is about 4.2 Km long with U-turns forming a loop. Trajectory #2 is about 2.4 Km long, and includes several roundabouts with U-turns forming a loop. These first two trajectories are shown respectively in blue and yellow in the left panel of Figure 6. Trajectory #3 is about 3.2 Km long, as shown in the right panel of Figure 6. The vehicle was driven along the selected trajectories for the learning phase and to perform the tests. The vehicle speed was between 0 and 60 km/h with a mean value around 30 km/h.

Regarding the environment, the Rambouillet area comprises tree canopies and surrounding vegetation. Dense trees and vegetation are present on many parts of trajectories #2 and #3. Buildings are present on some parts of trajectory #1 while other parts can be considered as open sky. For the purposes of this paper, the three trajectories will be considered as belonging to the same environment (sub-urban), since a finer classification is difficult. The HD map used has centimeter accuracy, but some markings are missing and others rendered obsolete by environment change. Regarding the quality of the lane markings, some have been painted recently, while others are older, but this is something for which we do not possess a usable classification. However, the onboard camera provides an indication of the quality of the detected lane marking. When the quality is below a threshold value, lane-marking information is discarded before the sensor fusion with FDE. Below, the reporting of our experiments is separated into two subsections, dealing respectively with internal integrity and external integrity. The first of the two subsections is a validation of the multi-sensor data fusion with the FDE step. The second is the validation of the PL computation based on Student's t-distribution. The strategy of learning the dof in the AT and CT directions is presented in detail.

A. Internal integrity

We first show experimentally the importance of having an FDE scheme in the data fusion process. Table I shows the errors before and after the FDE step for the three test trajectories. A reduction can be seen in the mean and max errors following the FDE step. For example, for the test trajectory #1, in the AT direction, the largest error falls from 20.42 m before the FDE step to 1.06 m after the FDE step. Likewise, in the CT direction, the error falls from 19.70 m to 1.03 m. The results for this trajectory are among the best results in our dataset in terms of accuracy as a result of the small GNSS errors in this experiment. Please note that for the test trajectory #2, the FDE does not significantly improve the accuracy because there were few outliers. For the trajectory #3, the improvement is mainly in the CT direction, where the mean error falls from 0.53 m before the FDE step to 0.46 m after the FDE step. The errors are in particular map errors detected at the beginning of the trajectory.

For a clearer illustration, the detailed results corresponding to the test trajectory #1 are shown below.

The residual used for the fault detection with the associated threshold is shown in figure 7. When this residual indicates the presence of faults, the set of residuals used for the fault exclusion is generated as shown in figure 8. Figure 8a shows the residual corresponding to the GNSS measurements, and figures 8b and 8c show the residuals corresponding to the camera measurements on the left and right sides, respectively. By comparing these residuals to the threshold value, it can be seen that the errors are mostly detected on the camera measurements. The differentiation between camera faults (including mapmatching) and map faults is shown in figure 9. Blue indicates that the error comes from the map, while red indicates that the origin of the error cannot be specified, for one of the two reasons presented in section V-A. Looking, for example, at figure 8b, it will be remarked that around the sampling instant 9693 (highlighted zone), r cameraL2 indicates the presence of faults while r cameraL1 is error free. Therefore, at these sampling instants, the faults are deemed to come from the map, and this is represented in figure 9a by the blue lines around the sampling instant 9693. In the same way, around the sampling instant 1522 (highlighted zone), the map is deemed erroneous, since r cameraL1 indicates the presence of faults while r cameraL2 is error free. Figure 10 shows two areas where a map fault has been detected. In these areas, the camera visualizes two lane markings on the left sides (the link border and the centerline) while only the link border is present in the map. This discrepancy led to an incorrect association. Therefore, the detected map faults are actually map errors.

Figure 11 shows the errors when using only GNSS measurements in relation to the multi-sensor data fusion before and after the FDE step. Before the FDE step, the effect of the errors on the position estimates is noticeable. The surges visible before the FDE step (in green) disappear after the FDE step (in red). The camera leads to a noticeable improvement in the CT direction with respect to the GNSS-only case. In the AT direction, the data fusion (with FDE) eliminates some outliers that occur with GNSS only, but does not improve overall accuracy, since the visual features do not provide longitudinal constraints. In some parts of the trajectory, the longitudinal error after the data fusion and after the FDE step is greater than with GNSS only, and this is mainly due to the drift of the odometry model. In other cases, some of the increases (the small surges) are due to undetected map-matching errors. This is an issue that could have been handled by optimizing the threshold, but that was not the purpose of this paper.

B. External integrity

In order to study the computation of the PLs, a TIR of 10 -3 was chosen, since we did not have enough samples to test with smaller integrity risks.

The procedure for tuning the best dof value is illustrated in table II. After the FDE step, the IR was computed on the training trajectories for different values of the dof (ν = 4 to ν = 100, corresponding to a Gaussian distribution). In the CT direction, the largest value of the dof that satisfies the IR requirement overall is seen to be ν = 9 (since 0.87×10 -3 <10 -3 ), although the TIR is not satisfied separately for the training trajectory #3.

The same methodology is applied in the AT direction, where ν = 5 has to be chosen. It will be remarked that this dof is not the same as in the CT direction, given that GNSS measurements are the only exteroceptive information available in the AT direction (the camera only detects marking in the CT direction). This explains the large error in the AT direction in some parts of the trajectory, where the GNSS is sometimes biased by about 3 m. This result would be different with a differential GNSS receiver, or if the GNSS biases were estimated in the filter. Moreover, it is clear from table II that it is unsuitable to use a Gaussian distribution for bounding the errors in this environment, since for ν = 100 the measured IRs are a long way above the specified values. To conclude, for this suburban environment and using the considered localization system, Student's tdistributions with ν = 9 and ν = 5 are required to bound the errors in the CT and AT directions respectively.

Figure 12 shows the errors with the computed PLs in the AT and CT directions after the FDE step applied to the test trajectory #1. The dof of the Student's t-distribution is set to ν = 5 in the AT direction and ν = 9 in the CT direction. Using these values, the computation of the PLs is consistent since the empirical IR is equal to 0 (which is smaller than 10 -3 ) in both directions. When no FDE step is added, and using the same dof as before, the measured frequencies of errors that exceed the PL (empirical IR) were 0.027 and 0.051 in the AT and CT directions respectively. Likewise, figure 13 shows similar results and highlights the importance of the FDE step for external integrity. Indeed, without FDE, the inconsistency of the computed PL persists in presence of large biases even when small degrees of freedom are used (for ν = 2.4, the empirical IRs were 0.019 and 0.009 in the AT and CT directions respectively). In addition, it will be noted that the PL decreases in the presence of undetected faults, since the erroneous measurement leads to an increase in the quantity of information, which lowers the values in the covariance matrix. Therefore, computing consistent PL in the presence of faulty measurements is not possible with the presented methodology.

Let us now look at the consistency of the computed PL on the different test trajectories (see Table III). In the CT direction (resp. AT direction), the Student's t-distribution with As an illustration of the values of the calculated PLs with the chosen dof, figure 14 shows their values on the test trajectory #2. Several rounds on the roundabouts explain the sinusoidal shape that appears in this figure.The camera cannot operate on highly curved roads, which explains the unavailability of its measurements on the roundabouts and the increase in the PL value, which is around 0.97 m in the CT direction when camera measurements are available, and about 3.4 m in the AT direction.

Figure 15a reports other results obtained with the test trajectory #3 with map errors. The intervals where PLs in the CT direction show an increase correspond to roundabouts. In this experiment, some map errors were automatically detected by the system. Figure 15b shows the residuals corresponding to the left camera measurements. From sampling instants 931 to 2278, the residual r cameraL2 indicates that the measurements C 0,left,2 is erroneous, while r cameraL1 indicates that C 0,left,1 is error-free. As discussed in section V, the error comes from the map. This was verified afterwards by comparing manually the environment and the map: the camera detected two markings on the left side, while only one features in the HD map. Despite map errors, the filter provides PLs that are always reliable. Finally, it should be noted that using a Gaussian distribution for the data fusion is a good approach, since the Gaussian distribution models the central part well. Since the test statistics are created according to the filter state estimates (considered as Gaussian), we retained the Chi-squared distribution for threshold setting. For the FDE, the false alarm probability is set to 0.05, which is large enough to be handled well by a Gaussian distribution. However, this would not be the case if the false alarm probability was set to 10 -3 .

VII. CONCLUSION AND FUTURE WORK

In this work, a general localization method for utilizing multi-lane marking detection with integrity considerations was proposed. It was tested using a GNSS receiver, DR measurements, a camera and an HD map. An FDE strategy based a bank of information filters was used to determine the origin of the faults thanks to a set of residual tests. We then proposed an accompanying method for bounding the errors, by replacing the Gaussian distribution with a Student's t-distribution. Our experimental results clearly show that this approach for PL computation is more relevant than classical methods that use Gaussian distributions. The protection levels are computed in the AT and CT directions in order to give a concrete meaning of PL for classical intelligent vehicle navigation tasks.

The performance of the approach was tested on different trajectories. The results obtained demonstrate the importance of the FDE stage and the characterization of the positioning errors by a heavy-tailed law for bounding errors when some faulty measurements remain in the computation process after the FDE step. Under the conditions of our experimental trials, we concluded that a Student's t-distribution with 9 dof is required in the CT direction, while 5 is required in the AT direction. We believe that the proposed method can be used in different land navigation environments. To adjust the dof of the t-distribution, data acquisitions with a ground truth system are necessary in the target environment (urban, suburban, open sky, etc.). This is in line with the strategy for deploying level 4 SAE autonomous vehicles within a well-defined operational design domain.

Figure 1 .

 1 Figure1. General architecture applied to the problem addressed in this paper. The fusion process takes as input information from GNSS, odometer, camera and HD map. The system performs an FDE step before providing pose estimates X k/k along with covariance matrices P k/k and PLs.

Figure 3 .

 3 Figure 3. Multi-sensor fusion with FDE. The filters use the camera measurements and the HD map as input.

Figure 4 .

 4 Figure 4. Comparison of two Student's t-distribution choices (same degree of freedom): in red: Gaussian distribution and the associated PL 1 , in blue: Student's t-distribution with shape matrix equal to the covariance of the Gaussian and the associated PL 2 , in green: Student's tdistribution with covariance matrix equal to the covariance of the Gaussian and the associated PL 3 .

Figure 5 .

 5 Figure 5. Display of the available HD map (in white) in the city of Rambouillet, France.

Figure 6 .

 6 Figure 6. The three different trajectories. The red color represents the backbone of the map. Left plot: In yellow, trajectory #1 (~4.2 km) with a round trip. In blue, trajectory #2 (~2.4 km) with a round trip and a number of roundabouts. Right plot: In yellow, trajectory #3 (~3.2 km). On the right side of this trajectory, the lane markings are not stored in the HD map but are detected by the camera.

Figure 7 .

 7 Figure 7. Residual r for fault detection. In green the threshold value (test trajectory #1).

Figure 8 .Figure 9 .

 89 Figure 8. Residuals r i for fault exclusion. In green the thresholds (test trajectory #1).

Figure 10 .

 10 Figure 10. In blue the areas where map errors are detected, and in yellow the missing markings.

Figure 11 .

 11 Figure 11. Errors after data fusion with FDE (in red) compared to data fusion without FDE (in green) and to GNSS (in blue). Sampling period of 0.02s corresponding to the period of the DR sensors. Test trajectory #1.

Figure 12 .

 12 Figure 12. PLs after the FDE step in the AT and CT directions. Test trajectory #1.Table III TEST STAGE. MEASURED IR APPLIED TO THE THREE TEST TRAJECTORIES. Measured IR Test Trajectories AT (ν = 5)

Figure 13 .

 13 Figure 13. PLs without FDE in the AT and CT directions for two Student's t-distributions (ν = 5 and ν = 2.4). Test trajectory #1.

Figure 14 .

 14 Figure 14. Test trajectory #2 :PLs in the AT and CT directions.

Figure 15 .

 15 Figure 15. Test trajectory #3: (a) PLs in the AT and CT directions, (b) left side camera residual with map errors visible at the beginning.

  global ENU frame R O and is attached to the body frame R B situated at the middle of the rear wheel axis. GNSS measurements are given in the antenna frame R G . The camera measurements are expressed in the frame R M situated at a distance P x in front of R B . The camera measures the lateral distance C 0 to the point L in a road marking segment [AB].

	𝐲 𝐎				𝒙 𝑴	𝑩
	𝑃 𝑥				𝐂 𝟎,𝐫𝐢𝐠𝐡𝐭,𝟏
				𝐑 𝐌
	𝑦 𝐵	𝑦 𝐺	𝑥 𝐵 𝐑 𝐆	𝑥 𝐺	𝒚 𝑴 𝑳
	𝒚	𝐑 𝐁 𝜃
			𝑨
	𝐑 𝐎
	𝒙			𝐱 𝐎
	Figure 2. Frames representation. The vehicle pose	x y θ	T is
	expressed in a

Table I HORIZONTAL

 I POSITION ERROR IN METERS BEFORE AND AFTER THE FDE STEP APPLIED TO THE TEST TRAJECTORIES.

			After FDE			Before FDE
	Test trajectories	error	max abs(AT)	max abs(CT)	error	max abs(AT)	max abs(CT)
	Traj. #1 (~4.2 km)	0.45	1.06	1.03	0.80	20.42	19.7
	Traj. #2 (~2.4 km)	2.48	3.28	2.93	2.48	3.29	2.93
	Traj. #3 (~3.2 km)	2.11	3.02	2.97	2.11	3.42	2.95

http://www.gnss-escape.eu/

http://www.sappart.net/

ACKNOWLEDGMENT

This work was carried out within SIVALab, joint laboratory between Renault and Heudiasyc UMR UTC/CNRS.