# Microscopic origin of magnetism and magnetic interactions in ferropnictides

M.D. Johannes, I.I. Mazin  Code 6393, Naval Research Laboratory, Washington,
D.C. 20375

(Printed on )

###### Abstract

One year after their initial discovery, two schools of thought have
crystallized regarding the electronic structure and magnetic properties of
ferropnictide systems. One postulates that these are itinerant weakly
correlated metallic systems that become magnetic by virtue of spin-Peierls
type transition due to near-nesting between the hole and the electron Fermi
surface pockets. The other argues these materials are strongly or at least
moderately correlated, the electrons are considerably localized and close to a
Mott-Hubbard transition, with the local magnetic moments interacting via
short-range superexchange. In this paper we argue that neither picture is
fully correct. The systems are moderately correlated, but with correlations
driven by Hund’s rule coupling rather than by the on-site Hubbard repulsion.
The iron moments are largely local, driven by Hund’s intra-atomic exchange.
Superexchange is not operative and the interactions between the Fe moments are
considerably long-range and driven mostly by one-electron energies of all
occupied states.

###### pacs:

Pacs Numbers:

Ferropnictides are still attracting widespread attention from researchers both
inside and outside the field of superconductivity. There is now a nearly
universal agreement that magnetism and, specifically, proximity to an
antiferromagnetic “stripe-order” transition plays a major role in the physics
of these compounds. There is also growing evidence that the magnetic
properties and correlation effects in this system are not controlled by the
Hubbard $U$ as in cuprates (spectroscopy tells us that Hubbard correlations
are weak and the effective $U$ is on the order of 1 eVPES , smaller than the
bandwidth; first principles calculations of $U$ support thisU ). On the other
hand, the multiband character of Fe bands and the large intra-atomic (Hund’s)
exchange coupling in Fe suggest that the Hund’s $J$ may play the main role in
the magnetism..

As opposed to the Hubbard $U,$ the Hund’s $J$ is generally well accounted for
in density-functional calculations (where it is called the Stoner $I$).)
Indeed, the local density approximation (LDA) correctly predicts the
particular antiferromagnetic and structural ground state of undoped
ferropnictides, in striking contrast to the cuprates. In view of this, it is
instrumental to trace down the origin of magnetism $within$ LDA, and to
disentangle the nature of magnetic interactions captured by this approach. It
is highly likely that the physics uncovered by density functional theory will
reflect the actual physics of these systems. Given the heated (but largely
devoid of solid facts) discussion of whether antiferromagnetism in pnictides
is due to Fermi surface nesting or to second neighbor superexchange (See Ref.
MS  for a review), a clear understanding of, at least, the message that LDA
calculations send seems highly necessary.

In this paper, we analyze the magnetic interactions and demonstrate that
neither of the above two views (often presented as an axiomatic dilemma) is
correct. The magnetism appears due to $local$ Hund’s rule coupling, while the
particular ground state is selected by itinerant, essentially one-electron
interactions, of which only a small part is played by the Fermi surface
nesting. We will also explain why conventional “Anderson-Kanamori”
superexchange is not operative here, and will show some striking examples
where calculations and experiment contradict both local superexchange and
spin-Peierls pictures, yet are perfectly understandable on the basis of one-
electron energy balance.

We start with a qualitative analysis. The Hund’s rule coupling energy in
density functional theory is expressed as $E_{H}=-Im^{2}/4,$ where $m$ is the
magnetization of an Fe ion, and $I\approx 0.8$ eV is the Stoner factor for Fe
(0.9 eV in GGA). Depending on the material, the self-consistent magnetic
moment on Fe appears to be between 1.5-2 $\mu_{B}$ in LDA and 1.8-2.5 in GGA.
The corresponding energy gain even in LDA is 0.5 eV, which is remarkably
large. In other words, every individual Fe wants to be strongly magnetic and
the advantage of spin polarization should lead to a magnetic ground state at
the mean field level, unless an unusually large kinetic energy penalty exists.
However, this is exactly the case for the formation of a ferromagnetic
configuration. To create a magnetization $m$ on Fe, one needs to move
approximately 1.15 (to account for the relative share of Fe-d orbitals at the
Fermi level) spin-minority electrons into unoccupied spin-majority states,
incurring an energy loss of $\approx(1.15m)^{2}/N_{\uparrow}(E_{F}).$ The
density of states (DOS) per Fe, $N_{\uparrow}(E_{F}),$ varies between 1 and
1.5 eV${}^{-1},$ depending on the system, creating an energy loss for $m=1.5$
$\mu_{B}$ of 0.5-0.8 eV. This cost is about as large as the Hund’s rule energy
gain estimated above. This shows that the system is on the verge of a
ferromagnetic instability, but nothing more.

In low-DOS metals, magnetization without a large cost in kinetic energy is
possible if some type of antiferromagnetic arrangement is formed (cf. metal Cr
and Mn). For a broad band metal, this narrows the conductivity band, but as
long as the exchange splitting is smaller than the bandwidth, the cost is
small. Because in ferropnictides the calculated bandwidth is 5-6 eV and the
exchange splitting $mI$ is at most 2 eV, this mechanism should be very
favorable.

It is interesting to consider how the system determines which particular AFM
arrangement is most profitable from the point of view of the one-electron
energy (note that LDA calculations can be forced to converge to nearly any AFM
pattern, but not to a ferromagnetic state). If the resulting magnetization is
small, the answer is obvious: the second derivative of the total energy with
respect to magnetization is defined by the noninteracting susceptibility at
the AF wave vector Q, $\partial^{2}E/\partial
m^{2}=-\chi_{0}^{-1}(\mathbf{Q)}$ (with the small caveat that an actual spin
density wave is not a single harmonic, but includes all wave vectors Q+G,
where G is a reciprocal lattice vector). The imaginary part of $\chi_{0}$ is
directly related to Fermi surface nesting, being defined, in the constant
matrix elements approximation, as
$\sum_{ij}\int\delta(\varepsilon_{\mathbf{k}i})\delta(\varepsilon_{\mathbf{k+Q,}i})d\mathbf{k,}$
while the (actually relevant) real part collects information from all states
and may or may not have any relation to the nesting conditions (for a detailed
discussion see Ref. CDW ).

Geometrical nesting, as a property of the Fermi surface, becomes even more
disconnected from a real instability in the strongly nonlinear regime,
$m\gtrsim 1$ $\mu_{B}$, which is the case for ferropnictides. Monitoring the
evolution of the electronic bands with increasing spin polarizationFSM , one
observes that at $m\sim 1$ $\mu_{B}$ the resulting bands can in no way be
described as anticrossing downfolded nonmagnetic bands with partial gapping of
the Fermi surface. Rather, the entire Fe d band is fully restructured.
Although the lowest-energy AFM state wave vector indeed coincides with the
quasi-nesting wave vector in some cases, it is not always true, as exemplified
by the case of FeTe that we discuss later.

It should be noted that while quasi-nesting is not particularly relevant for
the long-range ordering in the undoped crystals, it does define the low-energy
excitations in non-magnetic phases and these can perfectly well mediate
superconductivity.

Figure 1: (color online) Top-down view of the a) checkerboard b) stripe and c)
doublestripe magnetic patterns for a single FeAs or FeTe layer. The light
colored sites have majority up spin and the darker sites have majority down
spin.

Having established a general framework, we now address specific examples.
First, we investigate checkerboard, stripe, and double-stripe magnetic
structures (See Fig. 1) and show that the stripe order is lower in energy than
either the checkerboard or the double-stripe structure for the 122 systems,
but not for FeTe. We use BaFe2As2 as an example, but the results for LaFeAsO
are very similar. Our calculations were performed using an all-electron, full-
potential LAPW package WIEN2k, in the Generalized Gradient Approximation,
similar to Ref. PRB . All structures were fully relaxed (except where stated
otherwise) using the Vienna Ab-Initio Simulation Program (VASP) vasp , with
the PAW formulation PAW  and also using GGA. In Table 1 we show the magnetic
stabilization energies of the three different antiferromagnetic structures.

Table 1: Stabilization energies for various magnetic configurations in the 122 and FeTe systems. All energies are per Fe atom. | checkerboard | stripe | double stripe
---|---|---|---
BaFe2As2 | 16 meV | 94 meV | 0.6 meV
FeTe | – | 207 meV | 230 meV

In Fig. 2a,b,c, we show the DOSs for BaFe2As2 in each of the three magnetic
configurations along with the nonmagnetic DOS. Compared to the nonmagnetic
DOS, we see that the checkerboard pattern has a very similar spectrum at and
near the Fermi energy and gains one-electron energy by shifting spectral
weight from the region between -0.5 and -1.0 downward to the region between
-1.0 and -2.0. The doublestripe pattern actually incurs an energy penalty at
and just below EF, but gains energy by shifting weight downward from between
-0.2 and -0.7 to between -1.0 and -2.0. The ground state configuration, in
contrast to the other two, gains energy all the way from EF to -0.9 by
shifting weight downward. This is accomplished through the opening of a large
pseudogap (this terminology has no connection with the pseudogap in cuprates
and simply signifies a depression in one-electron DOS around the Fermi level).
Though all three magnetic configurations are stable with respect to a
nonmagnetic state, it is visibly the case that the stripe ordering has the
greatest one electron energy advantage. This is reflected in the much larger
gain in total energy (See Table 1).

Figure 2: (color online) The densities of states for BaFe2As2 in the non-
magnetic configuration in comparison to a) checkerboard magnetic pattern b)
stripe (ground state) magnetic pattern and c) double stripe magnetic pattern

Let us now compare the results with the same calculations for FeTe. As
indicated in a number of papers, FeTe is always formed with an excess Fe, so
the fact that experiment gives the double stripe structure as the low-
temperature ground state Li  should be taken $cum$ $grano$ $salis.$ However,
as Table 1 shows, it is definitely the stoichiometric ground state in density
functional calculations, and this is the only thing that matters for our
analysis note1 . We note here that we do not fully relax the FeTe structure,
but only relax the internal positions. As before these relaxations are done
separately for magnetic and nonmagnetic cases.

Figure 3: (color online) The densities of states for FeTe in the non-magnetic
configuration in comparison to a) stripe magnetic pattern and b) double stripe
magnetic pattern.

Table 1 indicates that for FeTe, as opposed to BaFe2As${}_{2},$ the energies
of the single and double stripe phases are relatively close. This suggests
that the crude method of determining the ground state by looking at the DOS
may not work here, as the DOSs for the two AFM phases will probably be
similar. Indeed, this is what we see in Fig. 3a,b where a large downward shift
of spectral weight is visible for both patterns. Interestingly, the
nonmagnetic Fermi surfaces in FeTe are extremely similar to those in the 122
and 1111 materials, whereas nonmagnetic DOS and the magnetic electronic
structure are quite different. This reinforces that FS nesting, which would be
nearly identical for BaFe2As2 and FeTe, is not driving the magnetic order.

Our calculations also provide a strong argument against superexchange. Looking
at the patterns in Fig. 1, it is easy to see that for both the stripe and
double stripe patterns, the first neighbor exchange, $J_{1}$, does not
contribute, due to equal numbers of aligned and anti-aligned spins. The second
neighbor exchange, $J_{2}$, would have to be stronger than $J_{1}/2$ in order
for the stripe pattern to be energetically favorable over the checkerboard
pattern in a superexchange picture. It has often been argued that this
situation is not unreasonable since the Fe-As-Fe paths available for $J_{1}$
and $J_{2}$ are similar. For double stripe order in the FeTe system, however,
both $J_{1}$ and $J_{2}$ cancel, leaving only $J_{3}$ to establish the
ordering. Considering the remarkably strong stability (compared to non-
magnetic) calculated for double stripe order (See Table 1), this is hard to
rationalize. Furthermore, the energy term for stripe is $J_{2}$ \- $J_{3}$
(compared to $J_{3}$ alone for double stripe). For double stripe to stabilize,
$J_{3}$ could be no smaller than $J_{2}/2$, but the ”similar hopping paths”
justification used for $J_{1}$ and $J_{2}$ and is not available: the third
neighbor exchange path is more than twice as long as the second neighbor one
and involves As-As hopping across a distance of a full lattice constant. Thus,
the existence and stability of the double stripe order severely strains the
credibility of the superexchange picture. This is, in fact, to be expected
since superexchange is not efficient when the bandwidth is much larger than
the energy cost of flipping an electron’s spin, which is precisely the case
here.

This does not, however, mean that one cannot map the dependence of the total
energy onto a suitable short-range exchange model. In fact, it is hard to
imagine a case in which this would not be possible. Yet, in carrying out this
procedure for ferropnictide systems, one should be aware of the following
caveats:

(1) There is no microscopic justification (as for instance in the Hubbard
model) for introducing any $J-t$ (or $J_{1}-J_{2}-t)$ Hamiltonian.

(2) There is no guarantee that this kind of mapping can be stopped at first or
second neighbors. In fact, accurate calculations show that at least some of
the exchange parameters in these mappings decay as $1/R^{3},$ just as in metal
ironY-A .

(3) The resulting exchange parameters strongly and qualitatively depend on the
long-range order established in the system. In particular, the parameters that
can be used to describe the ordered state cannot be used to describe the spin
fluctuations, and vice versa. (See Ref. Y-A  and references therein.)

(4) In the absence of superexchange, there is no reason to believe that the
total energy can be mapped onto a Heisenberg model,
$\sum_{ij}J_{ij}\mathbf{S}_{i}\cdot\mathbf{S}_{j}.$ In fact, direct
calculations show that at least one biquadratic term needs to be added to map
the total energy onto the mean-field Hamiltonian, $\sum
K(\mathbf{S}_{i}\cdot\mathbf{S}_{i+\mathbf{1}})^{2},$ where $K\sim J$. Ole

Figure 4: The Fermi surfaces of stripe-ordered BaFe2As2. Top panel shows the
’reverse distortion’ in which the Fe-Fe distance is lengthened along the like
spin direction and shortened along the unlike spin direction. The bottom
panels shows the fully relaxed calculation which reproduces the experimentally
observed distortion (to within a few percent).

We now switch our attention to the structural transition observed
simultaneoulsy with the magnetic one in the 122 systems. Density functional
calculations very accurately reproduce the experimentally observed distortion
in which Fe ions along the stripe direction are closer to one another than Fe
atoms belong to adjacent stripes yildirim ; jesche . We investigated whether
the structural distortion, like the magnetic ordering, can be understood in
terms of one electron energies by calculating the DOS for a variety of small
changes in the $a$ and $b$ lattice constants. In contrast to changing the
magnetic pattern, changing the structural distortion has very little effect on
the DOS away from the Fermi energy. There were no large shifts of spectral
weight to lower energies, though small shifts of the order of 0.05 eV did
occur and these were within 0.5 eV of the Fermi energy (for comparison see the
heavy restructuring of the DOS in Figs. 2 and 3). The distortion can therefore
be treated as a linear perturbation with a one-electron energy lowering
observable at (or very near) the Fermi energy. Specifically, we find that the
lowest energy structure corresponds to the smallest Fermi surface area. As an
example, in Fig. 4 we show the Fermi surface in the magnetic Brillouin zone of
the fully relaxed (lowest energy) structure and a ’reverse distortion’ in
which the distances between like and unlike spins are reversed from the
correct configuration. The change in the size of the Fermi surface is clearly
visible. We were unable to engineer a further minimization of the Fermi
surface with any choice of in-plane distortions other than the optimal energy
one.

In conclusion, we have shown that the relevant physics with respect to the
magnetic ordering and structural distortion in the ferropnictides lies in the
one-electron energies. Our results resolve the superficially binary choice
between superexchange and Fermi surface nesting in favor of a third mechanism
that is neither fully localized nor fully itinerant. One-electron energy is
gained throughout an energy range of at least 1 eV below EF and the ground
state is determined by which magnetic pattern most effectively exploits a
downshift in spectral weight, not by fermiology. On the other hand, the Fermi
surface itself is the operative feature for determination of the structural
distortion. The energy minimum for an in-plane distortion corresponds to a
simultaneous minimization of the Fermi surface area.

We thank H. Eschrig, K. Koepernik and T. Yildirim for useful and engaging
discussions related to this work. We acknowledge funding from the Office of
Naval Research.

## References

  * (1) T. Kroll, S. Bonhommeau, T. Kachel, H.A. Duerr, J. Werner, G. Behr, A.Koitzsch, R. Huebel, S. Leger, R. Schoenfelder, A. Ariffin, R. Manzke, F.M.F. de Groot, J. Fink, H. Eschrig, B. Buechner, M. Knupfer, Phys. Rev. B 78, 220502 (2008); F. Bondino, E. Magnano, M. Malvestuto, F. Parmigiani, M. A. McGuire, A. S. Sefat, B. C. Sales, R. Jin, D. Mandrus, E. W. Plummer, D. J. Singh, N. Mannella; Phys. Rev. Lett. 101, 267001 (2008); E. Z. Kurmaev, R. G. Wilks, A. Moewes, N. A. Skorikov, Yu. A. Izyumov, L. D. Finkelstein, R. H. Li, X. H. Chen, Phys. Rev. B 78, 220503(R) (2008); V. I. Anisimov, E. Z. Kurmaev, A. Moewes, I. A. Izyumov, Physica C (Special Issue), to be published; T.D. Devereaux et al, to be published.
  * (2) V. I. Anisimov, Dm. M. Korotin, M. A. Korotin, A. V. Kozhevnikov, J. Kuneš, A. O. Shorikov, S. L. Skornyakov, S. V. Streltsov, J. Phys.: Condens. Matter 21, 075602 (2009); T. Miyake, L. Pourovskii, V. Vildosola, S. Biermann, A. Georges, J. Phys. Soc. Jpn. 77,  Supplement C, 99 (2008).
  * (3) I.I. Mazin and J. Schmalian, Physica C (Special Issue), to be published.
  * (4) M.D. Johannes and I.I. Mazin, Phys. Rev. B77, 16535 (2008).
  * (5) M.A. Korotin, S.V. Streltsov, A.O. Shorikov, and V.I. Anisimov, JETP 107, 649 (2008).
  * (6) K.D. Belashchenko, V.P. Antropov, Phys. Rev. B78, 212505 (2008); T. Yildirim, Physica C (Special Issue), to be published.
  * (7) A. N. Yaresko, G.-Q. Liu, V. N. Antonov, O.K. Andersen, arXiv:0810.4469.
  * (8) I. I. Mazin, M. D. Johannes, L. Boeri, K. Koepernik, D. J. Singh, Phys. Rev. B78, 085104 (2008).
  * (9) G. Kresse, J. Furthmuller, Phys. Rev. B 54, 169 (1996).
  * (10) P. E. Blochl, Phys. Rev. B 50, 953 (1994).
  * (11) S. Li, C. de la Cruz, Q. Huang, Y. Chen, J.W. Lynn, J. Hu, Y-L. Huang, F-C. Hsu, K-W. Yeh, M-K. Wu, P. Dai, Phys. Rev. B 79, 054503 (2009).
  * (12) We are interested in how the double stripe structure manifests in stoichiometric FeTe within LDA. We cannot make meaningful comparisons with experiment and we therefore do not fully relax the FeTe structure, but only relax the internal positions. As before, relaxations are done separately for magnetic and nonmagnetic cases. We also neglect the checkerboard case as being of only secondary interest for this system.
  * (13) T. Yildirim, arXiv:0805.2888 (2008).
  * (14) A. Jesche, N. Caroca-Canales, H. Rosner, H. Borrmann, A. Ormeci, D. Kasinathan, K. Kaneko, H. H. Klauss, H. Luetkens, R. Khasanov, A. Amato, A. Hoser, C. Krellner, C. Geibel, Phys. Rev. B 78, 180504R (2008)


# A Particle Filtering Framework for Integrity Risk of GNSS-Camera Sensor
Fusion

Adyasha Mohanty, Shubh Gupta and Grace Xingxin Gao
Stanford University

## BIOGRAPHIES

Adyasha Mohanty is a graduate student in the Department of Aeronautics and
Astronautics at Stanford University. She graduated with a B.S. in Aerospace
Engineering from Georgia Institute of Technology in 2019.

Shubh Gupta is a graduate student in the Department of Electrical Engineering
at Stanford University. He received his B.Tech degree in Electrical
Engineering with a minor in Computer Science from the Indian Institute of
Technology Kanpur in 2018.

Grace Xingxin Gao is an assistant professor in the Department of Aeronautics
and Astronautics at Stanford University. Before joining Stanford University,
she was faculty at University of Illinois at Urbana-Champaign. She obtained
her Ph.D. degree at Stanford University. Her research is on robust and secure
positioning, navigation and timing with applications to manned and unmanned
aerial vehicles, robotics, and power systems.

## ABSTRACT

Adopting a joint approach towards state estimation and integrity monitoring
results in unbiased integrity monitoring unlike traditional approaches. So
far, a joint approach was used in Particle RAIM [1] for GNSS measurements
only. In our work, we extend Particle RAIM to a GNSS-camera fused system for
joint state estimation and integrity monitoring. To account for vision faults,
we derive a probability distribution over position from camera images using
map-matching. We formulate a Kullback-Leibler Divergence [2] metric to assess
the consistency of GNSS and camera measurements and mitigate faults during
sensor fusion. The derived integrity risk upper bounds the probability of
Hazardously Misleading Information (HMI). Experimental validation on a real-
world dataset shows that our algorithm produces less than 11 m position error
and the integrity risk over bounds the probability of HMI with 0.11 failure
rate for an 8 m Alert Limit in an urban scenario.

##  1 INTRODUCTION

In urban environments, GNSS signals suffer from lack of continuous satellite
signal availability, non line-of-sight (NLOS) errors and multi-path effects.
Thus, it is important to quantify the integrity or measure of trust in the
correctness of the positioning solution provided by the navigation system.
Traditional integrity monitoring approaches [3] provide point positioning
estimates i.e. the state estimation algorithm is assumed to be correct and
then the integrity of the estimated position is assessed. However, addressing
state estimation and integrity monitoring separately does not capture the
uncertainty in the state estimation algorithm. As a result, the integrity
monitoring becomes biased by the acquired state estimate leading to subsequent
faulty state estimation.

Recently, an approach towards joint state estimation and integrity monitoring
for GNSS measurements was proposed in Particle RAIM [1]. Instead of producing
point positioning estimates, Particle RAIM uses a particle filter to form a
multi-modal probability distribution over position, represented as particles.
Traditional RAIM [4] is used to assess the correctness of different ranging
measurements and the particle weights are updated to form the distribution
over the position. From the resulting probability distribution, the integrity
risk is derived using an approximate upper bound to the probability of HMI or
the reference risk. By incorporating the correctness of different measurement
subsets directly into the state estimation, Particle RAIM is able to exclude
multiple faults in GNSS ranging measurements. However, due to large errors
from GNSS measurements, Particle RAIM requires employing conservative measures
such as large Alert Limits to adequately bound the reference risk.

For urban applications, improved positioning accuracy from Particle RAIM is
necessary to provide adequate integrity for smaller Alert Limits. Since
measurements from GNSS are not sufficient to provide the desired accuracy, it
is helpful to augment GNSS with additional sensors that increase redundancy in
measurements. Sensors such as cameras are effective complimentary sensors to
GNSS. In urban regions, cameras have access to rich environmental features [5]
[6] [7] and provide superior sensing than GNSS which suffers from multi-path
and NLOS errors [3] [8] [9] [10].

Thus, with added vision, we need a framework to provide integrity for the
fused GNSS-camera navigation system to account for two categories of faults.
The first category includes data association errors across images, where
repetitive features are found in multiple images creating ambiguity during
feature and image association. This ambiguity is further amplified due to
variations in lighting and environmental conditions. The second category
comprises errors that arise during sensor fusion of GNSS and camera
measurements. Ensuring that faults in either measurement do not dominate the
sensor fusion process is paramount for maximizing the complimentary
characteristics of GNSS and camera.

Many works provide integrity for GNSS-camera fused systems utilizing a Kalman
Filter [11] framework or an information filter [12]. Vision Aided-RAIM [13]
introduced landmarks as pseudo-satellites and integrated them into a linear
measurement model alongside GPS observations. In [14], the authors implemented
a sequential integrity monitoring approach to isolate single satellite faults.
The integrity monitor uses the innovation sequence output from a single Kalman
filter to derive a recursive expression of the worst case failure mode slopes
and to compute the protection levels (PL) in real-time. An Information Filter
(IF) is used in [15] for data fusion wherein faults are detected based on the
Kullback-Leibler divergence (KL divergence) [2] between the predicted and the
updated distributions. After all detected faulty measurements are removed, the
errors are modeled by a student’s t distribution to compute a PL. A student’s
t distribution is also used in [16] alongside informational sensor fusion for
fault detection and exclusion. The degree of the distribution is adapted in
real-time based on the computed residual from the information filter. A
distributed information filter is proposed in [17] to detect faults in GPS
measurement by checking the consistency through log-likelihood ratio of the
information innovation of each satellite. These approaches model measurement
fault distributions with a Gaussian distribution although for camera
measurements, the true distribution might be non-linear, multi-modal, and
arbitrary in nature. Using a simplified linear measurement probability
distribution renders these frameworks infeasible and unreliable for safety-
critical vision augmented GNSS applications.

Another line of work builds on Simultaneous Localization and Mapping (SLAM)
based factor graph optimization techniques. Bhamidipati et al [5] derived PL
by modeling GPS satellites as global landmarks and introducing image pixels
from a fish-eye camera as additional landmarks. The raw image is categorized
into sky and non-sky pixels to further distinguish between LOS and NLOS
satellites. The overall state is estimated using graph optimization along with
an M-estimator. Although this framework is able to exclude multiple faults in
GPS measurements, it is not extendable to measurements from forward or rear
facing cameras that do not capture sky regions. Along similar lines,
measurements from a stereo camera along with GNSS pseudoranges are jointly
optimized in a graph optimization framework in [18]. GNSS satellites are
considered as feature vision points and pose-graph SLAM is applied to achieve
a positioning solution. However, graph optimization approaches also share the
same limitation as Kalman Filter based approaches: They produce point
positioning estimates and do not account for the uncertainty in state
estimation that biases integrity monitoring.

Overall, existing integrity monitoring algorithms for GNSS- camera fusion have
the following limitations:

  * 1

They address state estimation and integrity monitoring separately, similar to
traditional RAIM approaches.

  * 2

They accommodate camera measurements within a linear or linearizable framework
such as KF, EKF, or IF and become infeasible when camera measurements are not
linearizable without loss of generality.

  * 3

There is no standard way in literature to quantify the uncertainty in camera
measurements directly from raw images.

  * 4

They use outlier rejection techniques to perform fault detection and exclusion
after obtaining the positioning solution. There is no framework that accounts
for faults both independently in GNSS and camera as well as the faults that
arise during sensor fusion.

In our work, we overcome the above limitations by proposing the following
contributions. This paper is based on our recent ION GNSS+ 2020 conference
paper [19].

  * 1

We jointly address state estimation and integrity monitoring for GNSS-camera
fusion with a particle filtering framework. We retain the advantages of
Particle RAIM while extending it to include camera measurements.

  * 2

We derive a probability distribution over position directly from images
leveraging image registration.

  * 3

We develop a metric based on KL divergence [2] to fuse probability
distributions obtained from GNSS and camera measurements. By minimizing the KL
divergence of the distribution from each camera measurement with respect to
the GNSS measurement distribution, we ensure that erroneous camera
measurements do not affect the overall probability distribution. Stated
otherwise, the divergence metric augments the shared belief over the position
from both sensor measurements by minimizing cross-contamination during sensor
fusion.

  * 4

We experimentally validate our framework on an urban environment dataset [20]
with faults in GNSS and camera measurements.

The rest of the paper is organized as follows. In Section 2, we describe the
overall particle filter framework for probabilistic sensor fusion. In Sections
3 and 4, we infer a distribution over position from GNSS and camera
measurements, respectively. Section 5 elaborates on the probabilistic sensor
fusion of GNSS and camera measurements along with the proposed KL divergence
metric. In Section 6, we describe the integrity risk bounding. Sections 7 and
8 shows the experimental setup and the results from experimental validation on
the urban environment dataset, respectively. In Section 9, we conclude our
work.

##  2 PARTICLE FILTER FRAMEWORK FOR PROBABILISTIC SENSOR FUSION

The distribution over the position inferred from GNSS and camera measurements
is multi-modal due to faults in a subset of measurements. Thus, to model such
distributions, we choose a particle filtering approach that further allows us
to keep track of multiple position hypotheses rather than a single position
estimate. Although a particle filtering approach was used in Particle RAIM
[1], the authors only considered GNSS ranging measurements. In our work, we
extend the framework to include measurements from a camera sensor. Figure 1
represents our overall framework. We add the camera and probabilistic sensor
fusion modules to the framework proposed in [1].

Figure 1: Particle filter framework with probabilistic sensor fusion of GNSS
and camera measurements and integrity risk bounding. The highlighted modules
represent our contributions.The GNSS and Risk Bounding Modules are adopted
from Particle RAIM [1].

Our framework consists of the following modules:

  * •

Perturbation and propagation: Using noisy inertial odometry from the IMU, we
generate a set of motion samples, each of which perturbs the previous particle
distribution in the propagation step.

  * •

GNSS module: This module from Particle RAIM [1] takes GNSS ranging
measurements from multiple satellites, some of which may be faulty and outputs
a probability distribution over position using a fault-tolerant weighting
scheme described in Section 3. The particles from the GNSS module are
propagated to the camera module to ensure that the distributions from GNSS and
camera share the same domain of candidate positions.

  * •

Camera module and synchronization with motion data: The camera module takes a
camera image and matches it to the images in a map database using image
registration to generate similarity scores. The underlying state of the best
matched image is extracted and propagated forward to the current GNSS time
stamp by interpolating with IMU odometry. This step ensures that the
probability distributions from camera and GNSS measurements are generated at
the same time stamps. Finally, we use a categorical distribution function to
transform the similarity scores into a probability distribution over position
hypotheses as described in Section 4.

  * •

Probabilistic sensor fusion: This module outputs a joint likelihood over
positions from GNSS and camera measurements after fusing them with the
proposed KL divergence metric in Section 5.1. Particles are resampled from the
current distribution with Sequential Importance Resampling [21].

  * •

Risk bounding: We adopt the risk bounding formulation proposed in [1] to
compute the integrity risk from the derived probability distribution over the
position domain. Using generalization bounds from statistical learning theory
[22], the derived risk bound is formally shown to over bound the reference
risk in Section 6.

We elaborate on the various modules of our framework in the following
sections.

##  3 GNSS MODULE- PARTICLE RAIM

A likelihood model for the GNSS measurements is derived using the mixture
weighting method proposed in Particle RAIM [1]. Instead of assuming
correctness of all GNSS measurements, the likelihood is modeled as a mixture
of Gaussians to account for faults in some measurements. Individual
measurement likelihoods are modeled as Gaussians with the expected
pseudoranges as means and variance based on Dilution of Precision(DoP). The
GMM [23] [24] is expressed as:

$L_{t}(m^{t})=\sum_{k=0}^{R}\gamma_{k}\mathcal{N}(m_{k}^{t}|\mu_{X}^{t,k},\sigma_{X}^{t,k});\sum_{k=0}^{R}\gamma_{k}=1,$
(1)

where $L_{t}(m^{t})$ denotes the likelihood of measurement $m$ at time $t$.
$\gamma$ denotes the measurement responsibility or the weights of the
individual measurement components and $R$ refers to the total number of GNSS
ranging measurements. $\mu$ and $\sigma$ represent the mean and the standard
deviation of each Gaussian component inferred from DOP. $X$ refers to the
collection of position hypotheses denoted by particles and $k$ is the index of
the number of Gaussians in the mixture. The weights are inferred with a single
step of the Expectation-Maximization (EM) scheme [25] as shown in Figure 2.

Figure 2: Two steps of the EM scheme used to derive the weight of each
Gaussian likelihood in the GMM. In the expectation step, the local vote for
each particle is computed based on the squared-normal voting on the normalized
residual for a particle obtained with traditional RAIM. The overall confidence
is inferred by normalizing the votes and pooling them using Bayesian maximum a
posteriori (MAP) estimation.

To avoid numerical errors due to finite precision, the log likelihood of the
likelihood model is implemented by extending the input space to include
additional copies of the state space variable, one for each GNSS measurement
[26]. The new likelihood is written as:

$P\left(m^{t}\middle|X^{t},\chi=k\right)=\gamma_{k}\mathcal{N}\left(m_{k}^{t}\middle|\mu_{x}^{t,k},\sigma_{x}^{t,k}\right)\
;\sum_{k=1}^{R}\gamma_{k}=1,$ (2)

where $\chi$ is an index that denotes the associated GNSS measurement with the
particle replica.

##  4 CAMERA MODULE

To quantify the uncertainty from camera images, we use a map-matching
algorithm that matches a camera image directly to an image present in a map
database. Our method is implemented in OpenCV [27] and comprises three steps
shown in Figure 3.

Figure 3:  Generating probability distribution over position from camera
images.

Each block is elaborated below.

  * •

Database Creation: We assume prior knowledge of the geographical region where
we are navigating. Based on GPS coordinates, we select images from the known
area using Google Street View Imagery. These images along with their
associated coordinates form the database. Features are extracted from these
images and stored in a key point-descriptor format.

  * •

Image Registration: After receiving a camera test image, we extract features
and descriptors with the ORB [28] algorithm. Although we experimented with
other feature extraction methods such as SIFT [29], SURF [30], and AKAZE [31],
ORB was found most effective for extracting descriptors from highly blurred
images. The descriptor vectors are clustered with a k-means algorithm [32] to
form a vocabulary tree [33]. Each node in the tree corresponds to an inverted
file, i.e., a file containing the ID-numbers of images in which a particular
node is found and the relevance of each feature to that image. The database is
then scored hierarchically based on Term Frequency Inverse Document Frequency
(TF-IDF) scoring [33], which quantifies the relevance of the images in the
database to the camera image. We refer to these scores as the similarity
scores. The image with the highest score is chosen as the best match and the
underlying state is extracted.

  * •

Probability generation after synchronization: After extracting the state from
the best camera image in the database, we propagate the state to the same time
stamp as the GNSS measurement. The raw vehicle odometry is first synchronized
with GNSS measurements using the algorithm in [20]. Using the time difference
between the previous and current GNSS measurements, we linearly interpolate
the extracted state with IMU motion data as shown below.

$x^{t}=x^{t-1}+v^{t-1}dt+0.5a^{t-1}\ dt^{2}$ (3)

where $x^{t}$ refers to the 3D position at epoch $t$, $dt$ refers to the time
difference between successive camera measurements, and $v$ and $a$ are the
interpolated IMU velocity and accelerations at epoch $t$.

Next, we compute the Euclidean distance between the interpolated state and the
current particle distribution from GNSS measurements to obtain new similarity
scores. This step ensures that the probability distributions computed from
camera and GNSS measurements share the same domain of candidate positions. A
SoftMax function takes the scores and outputs a probability distribution over
position. Normalization of the scores enforces a unit integral for the
distribution.

$Q(n^{t}|X^{t})=\frac{\exp(\omega^{t})}{\sum_{c}\exp(\omega_{c}^{t})}$ (4)

where $Q$ is the probability distribution associated with camera measurement
$n$ at time $t$ over the position domain $X$, $\omega_{c}^{t}$ represents
computed distance score, and $c$ is the index for individual particles.

##  5 PROBABILISTIC SENSOR FUSION

After obtaining the probability distributions from GNSS and camera, we need to
form a joint distribution over the position. However, we need to ensure that
faults in camera measurements do not degrade the distribution from GNSS
measurements, one that is coarse but correct since the distribution accounts
for faults in the ranging measurements through the RAIM voting scheme. Thus,
we need a metric to identify and exclude faulty camera measurements leveraging
knowledge of the distribution from GNSS. Additionally, the metric should
assess the consistency of the probability distribution from each camera
measurement with respect to the GNSS distribution and mitigate inconsistent
distributions that result from vision faults. The KL divergence [34]
represents one way to assess the consistency of two probability distributions.
By minimizing the divergence between the distributions inferred from camera
and GNSS, we ensure that both distributions are consistent.

###  5.1 Kl Divergence: Metric Formulation

We provide a background on KL divergence prior to explaining our metric.

The KL divergence [34] between two discrete probability distributions, $p$ and
$q$, in the same domain is defined as:

$D_{KL}(p||q)=\sum\nolimits_{z\in\zeta}p_{z}\ log\ \frac{p_{z}}{q_{z}}$ (5)

where $\zeta$ represents the domain of both distributions and $z$ is each
element of the domain. In our work, we ensure that distributions from GNSS and
camera share the same position domain by propagating the particles from the
GNSS distribution to the camera module prior to generating the distribution
from camera measurements. Two important properties of the KL divergence are:

  * •

The KL divergence between two distributions is always non-negative and not
symmetrical [34]

$D_{KL}(p||q)\neq D_{KL}(q||p)$ (6)

where $D_{KL}(q||p)$ is the reverse KL divergence between the distributions
$p$ and $q$.

  * •

$D_{KL}(p||q)$ is convex in the pair $(p||q)$ if both distributions represent
probability mass functions (pmf) [34].

Leveraging the above properties, we formulate our metric below.

  * •

Mixture of Experts (MoE): We form a mixture distribution to represent
probability distributions from successive camera measurements, where a non-
Gaussian probability distribution is derived from a single camera image. Each
measurement is assigned a weight to represent its contribution in the mixture.
Instead of setting arbitrary weights, we leverage the GNSS distribution to
infer weights that directly correspond to whether a camera measurement is
correct or faulty. Thus, highly faulty camera measurements are automatically
assigned low weights in the MoE. The mixture distribution is given as:

$Q^{*}(n^{t}|X^{t})=\sum\limits_{j=1}^{K}\alpha_{j}^{*}\
Q^{j}(n_{j}^{t}|X^{t});\sum\limits_{j=1}^{K}\alpha_{j}^{*}=1$ (7)

where $Q^{*}(n^{t}|X^{t})$ represents the mixture distribution formed using
$K$ camera images between two successive GNSS time epochs.
$Q^{j}(n_{j}^{t}|X^{t})$ is the likelihood of a single camera image
$n_{j}^{t}$ recorded at time $t$ with $\alpha_{j}^{*}$ as the normalized
weight. $X^{t}$ are the particles representing position hypothesis and $j$ is
the index for the camera images. The weights are normalized below to ensure
that the MoE forms a valid probability distribution:

$\alpha_{j}^{*}=\frac{\alpha_{j}}{\sum\limits_{r=1}^{K}\alpha_{r}}$ (8)

where $\alpha_{j}^{*}$ is the normalized weight, $\alpha_{j}$ is the weight
prior to normalization, $r$ is the index for the number of camera images
between two successive GNSS time epochs, and $K$ is the total number of camera
measurements.

  * •

Setup KL divergence: We set up a divergence minimization metric between the
distributions from each camera measurement and all GNSS measurements.

${KL}_{j}\ ((\alpha_{j}\ Q^{j}\left(n_{j}^{t}\middle|X^{t}\right)\ ||\ P\
\left(m_{k}^{t}\middle|X^{t},\chi=k\right))=\sum_{i=1}^{S}{\left(\alpha_{j}\
Q^{j}(n_{j}^{t}|\ X^{t})\right)\ log\left[\frac{\left(\alpha_{j}\
Q^{j}(n_{j}^{t}|\ X^{t})\right)}{P\
\left(m_{k}^{t}\middle|X^{t},\chi=k\right)}\right]}$ (9)

where $||\ $ denotes the divergence between both probability distributions,
$S$ represents the total number of particles or position hypotheses across
both distributions, and $i$ is the index for the particles. $\ P\
\left(m_{k}^{t}\middle|X^{t},\chi=k\ \right)$ is the probability distribution
at epoch $t$ from GNSS measurements as defined in Equation (2), $\alpha_{j}\ $
is the unnormalized weight, and $j$ is the index for the camera measurement.

  * •

Minimize divergence: Using the convexity of the KL divergence (Property 2), we
minimize each divergence metric with respect to the unknown weight assigned to
the likelihood of each camera measurement. We abbreviate $\
P\left(m_{i}^{t}\middle|X^{t},\chi=i\ \right)$ as $P(x_{i})$ and
$Q\left(n_{j}^{t}\middle|X^{t}\right)\ $ as $Q(x_{i})$ for brevity and expand
Equation (9). Since $\alpha_{j}$ is independent of the summation index, we
keep it outside the summation and simplify our expansion below.

${KL}_{j}(Q|\left|P\right)=\ \alpha_{j}\sum_{i\ =\
1}^{S}Q\left(x_{i}\right)log\ \alpha_{j}\ +\ \alpha_{j}\sum_{i\ =\
1}^{S}Q\left(x_{i}\right)log\ Q\left(x_{i}\right)\ -\ \alpha_{j}\sum_{i\ =\
1}^{S}Q\left(x_{i}\right)log\ P\ \left(x_{i}\right)$ (10)

Taking the first derivative with respect to $\alpha_{j}$ we obtain,

${min}_{\alpha_{j}}{KL}_{j}\ \ (Q|\left|P\right)=\ log\ \alpha_{j}\sum_{i\ =\
1}^{S}Q\left(x_{i}\right)\ +\ \sum_{i\ =\ 1}^{S}Q\left(x_{i}\right)\ +\sum_{i\
=\ 1}^{S}Q\left(x_{i}\right)log\ Q\left(x_{i}\right)\ -\ \sum_{i\ =\
1}^{S}Q\left(x_{i}\right)log\ P\left(x_{i}\right)$ (11)

Equating the expression on the right to 0 and solving for $\alpha_{j}$ gives
us:

$\alpha_{j}=e^{k}\ ;\ k=\frac{\sum_{i=1}^{S}Q\left(x_{i}\right)\ log\
\frac{P\left(x_{i}\right)}{Q\left(x_{i}\right)}}{\sum_{i=1}^{S}Q\left(x_{i}\right)}-1$
(12)

We also perform a second derivative test to ensure that the $\alpha_{j}$ value
inferred is a minimum value of the divergence measure. Since the exponential
function with the natural base is always positive, $\alpha_{j}$ is always
positive as well. Thus, evaluating the second derivative gives us a positive
value.

$\frac{1}{\alpha_{j}}\sum_{i=1}^{S}Q(x_{i})>0$ (13)

  * •

Joint probability distribution over position: After obtaining the weights, we
normalize them using Equation (8). We obtain the joint distribution assuming
that the mixture distribution from camera measurements and the GMM from GNSS
measurements are mutually independent. The joint distribution is given as:

$P^{\ast}\left(n^{t},\ m^{t}\middle|X^{t}\right)=\
P\left(m_{i}^{t}\middle|X^{t},\chi=k\ \right)\
Q^{\ast}\left(n^{t}\middle|X^{t}\right)$ (14)

where $\ P\ \left(m_{k}^{t}\middle|X^{t},\chi=k\ \right)$ is the probability
distribution from GNSS measurements in Equation (2). We take the log
likelihood of the joint distribution to avoid finite precision errors.

##  6 INTEGRITY RISK BOUNDING

We upper bound the probability of HMI using the risk bounding framework
introduced in [1]. For a single epoch, the probability of HMI for a given
Alert Limit $r$ is defined as:

$R_{x*}(\pi)=\mathop{\mathbb{E}}_{x\sim\pi}[P(\|x-x^{*}\|\geq r)]$ (15)

where $R_{x*}(\pi)$ is the probability of HMI with reference position $x^{*}$
and mean distribution in position space induced by all posterior distributions
$\pi$. The distributions are created by generating samples around the measured
odometry and then perturbing the initial particle distribution. From the PAC-
Bayesian [35] formulation and as shown in [1], the reference risk
$\mathop{\mathbb{\textbf{R}}(\pi^{t})}$ upper bound is:

$\mathop{\mathbb{\textbf{R}}(\pi^{t})}\leq\mathop{\mathbb{\textbf{R}}_{M}(\pi^{t})}+\mathcal{D}_{Ber}^{-1}(\mathop{\mathbb{\textbf{R}}_{M}(\pi^{t})},\epsilon)$
(16)

The first and second terms refer to empirical and divergence risk,
respectively. We explain the computation of each term below.

The empirical risk $\mathop{\mathbb{\textbf{R}}_{M}(\pi^{t})}$ is computed
from a finite set of perturbed samples of size $M$.

$\mathop{\mathbb{\textbf{R}}_{M}(\pi^{t})}=\frac{1}{M}\sum_{i=1}^{M}\mathop{\mathbb{E}}_{x\sim\pi^{t}}[l(x,\pi_{u}^{t})],$
(17)

where, $l(x,\pi_{u}^{t})$ is the classification loss with respect to a motion
sample resulting in the posterior distribution being classified as hazardous.
$\pi$ refers to the mean posterior distribution at time $t$.

The divergence risk term
$\mathcal{D}_{Ber}^{-1}(\mathop{\mathbb{\textbf{R}}_{M}(\pi^{t})},\epsilon)$
accounts for uncertainty due to perturbations that are not sampled. First, we
compute the gap term $\epsilon$ using KL divergence [2] of the current
distribution from the prior and a confidence requirement in the bound
$\delta$.

$\epsilon=\frac{1}{M}(KL(\pi^{t}||\pi^{t-1})+log(\frac{M+1}{\delta}))$ (18)

where $\delta$ refers to the bound gap. The means of the prior and current
distributions are taken as $\pi^{t-1}$ and $\pi^{t}$. The prior and current
distributions are approximated as multivariate Gaussian distributions.

The Inverse Bernoulli Divergence [1] $\mathcal{D}_{Ber}^{-1}$ is defined as:

$\mathcal{D}_{Ber}^{-1}(q,\epsilon)=t\;\;s.t.\;\;\mathcal{D}_{Ber}(q||q+t)=\epsilon$
(19)

where $q||q+t$ is the KL divergence [2] between $q$ and $q+t$ and $q$ is given
by the empirical risk term. Finally, the Inverse Bernoulli Divergence [1] is
obtained approximately as:

$\mathcal{D}_{Ber}(q,\epsilon)=\sqrt{\frac{2\epsilon}{\frac{1}{q}+\frac{1}{1-q}}}$
(20)

##  7 EXPERIMENTS

###  7.1 Datasets

We test our framework on a 2.3 km long urban driving dataset from Frankfurt
[20]. We use GNSS pseudorange measurements, images from a forward-facing
camera, ground truth from a NovAtel receiver, and odometry from the IMU. The
dataset contains NLOS errors in GNSS measurements and vision faults due to
variations in illumination. In addition to the real-world dataset, we create
emulated datasets by inducing faults in GNSS and vision measurements with
various controlled parameters.

###  7.2 Experimental Setup and Parameters

  * •

Real-world dataset: We use GNSS ranging measurements with NLOS errors. For
simplicity, we estimate the shared clock bias by subtracting the average
residuals with respect to ground truth from all GNSS pseudoranges at one time
epoch.

  * •

Emulated dataset: First, we vary the number of satellites with NLOS errors by
adding back the residuals to randomly selected satellites. This induces clock
errors in some measurements which are perceived as faults. Secondly, we remove
the NLOS errors from all measurements but add Gaussian bias noise to
pseudorange measurements from random satellites at random time instances. The
number of faults are varied between 2-9 out of 12 available measurements at
any given time step. We induce faults in camera measurements by adding
blurring with a 21x21 Gaussian kernel and occlusions of 25-50 % height and
width to random images.

During the experimental simulation, a particle filter tracks the 3D position
(x,y,z) of the car and uses faulty GNSS and camera measurements along with
noisy odometry. Probability distributions are generated independently from
GNSS and camera and fused with the KL divergence metric to form the joint
distribution over positions. At each time epoch, the particle distribution
with the highest total log-likelihood is chosen as the estimated distribution
for that epoch. The integrity risk is computed from 10 posterior distributions
of the initial particle distribution and the reference risk is computed with
ground truth. Our experimental parameters are listed in Table 1.

Table 1: Experimental Parameters for Validation with Real-world and Emulated Datasets Parameter | Value | Parameter | Value
---|---|---|---
No. of GNSS measurements | 12 | Added Gaussian bias to GNSS measurements | 20- 200 m
No. of faults in GNSS measurements | 2-9 | No. of particles | 120
Measurement noise variance | 10 m 2 | Filter propagation variance | 3 m 2
Alert Limit | 8, 16 m | No. of odometry perturbations | 10

###  7.3 Baselines and Metrics

We use Particle RAIM as the baseline to evaluate our algorithm’s performance
for state estimation. The metric for state estimation is the root mean square
error (RMSE) of the estimated position with respect to ground truth for the
entire trajectory. The risk bounding performance is evaluated with metrics
derived from a failure event, i.e., when the derived risk bound fails to upper
bound the reference risk. The metrics are the following: failure ratio(the
fraction of cases where the derived risk bound fails to upper bound the
reference risk), failure error(the mean error during all failure events), and
the bound gap(average gap between the derived integrity risk) and the
reference risk.

For evaluating the integrity risk, we specify a performance requirement that
the position should lie within the Alert Limit with at least 90% probability.
A fault occurs if the positioning error exceeds the Alert Limit. The metrics
for integrity risk are reported based on when the system has insufficient
integrity or sufficient integrity [36], which respectively refer to the states
when a fault is declared or not. The false alarm rate equals the fraction of
the number of times the system declares insufficient integrity in the absence
of a fault. The missed identification rate is defined as the fraction of the
number of times the system declares sufficient integrity even though a fault
is present.

##  8 RESULTS

###  8.1 State Estimation

First, we test our algorithm with NLOS errors in GNSS ranging measurements and
added camera faults. Quantitative results in Table 2 demonstrate that our
algorithm produces 3D positioning estimates with overall RMSE of less than 11
m. Additionally, our algorithm reports lower errors compared to Particle RAIM
for all test cases. Our algorithm is able to compensate for the residual
errors from Particle RAIM by including camera measurements in the framework.
This leads to improved accuracy in the positioning solution.

Table 2: RMSE in 3D Position with NLOS errors and added vision faults No. of faults out of 12 available GNSS measurements | Particle RAIM-Baseline (meter) | Our Algorithm (meter)
---|---|---
2 | 18.1 | 6.3
4 | 19.1 | 6.1
6 | 16.9 | 5.9
9 | 26.6 | 10.6

For qualitative comparison, we overlay the trajectories from our algorithm on
ground truth and highlight regions with positioning error of greater than 10 m
in Figures 4 and 5. Trajectories from Particle RAIM show large deviations from
ground truth in certain regions, either due to poor satellite signal
availability or high NLOS errors in the faulty pseudorange measurements.
However, similar deviations are absent from the trajectories from our
algorithm which uses both GNSS and camera measurements. Our KL divergence
metric is able to mitigate the errors from vision and the errors from cross-
contamination during sensor fusion, allowing us to produce lower positioning
error.

(a) Particle RAIM (Baseline)

(b) Our Algorithm

Figure 4: State estimation under NLOS errors for 6 faulty GNSS pseudo range
measurements and added vision faults. Regions with positioning error greater
than 10 m are highlighted in red.

(a) Particle RAIM (Baseline)

(b) Our Algorithm

Figure 5: State estimation under NLOS errors for 9 faulty GNSS pseudo range
measurements and added vision faults. Regions with positioning error greater
than 10 m are highlighted in red.

Secondly, we test our algorithm with the emulated datasets. Quantitatively, we
plot the RMSE as a function of the added Gaussian bias value in Figure 6 and
as a function of the number of faulty GNSS ranging measurements in Figure 7.
For all validation cases, our algorithm produces an overall RMSE less than 10
m. Similar to the results from the real-world dataset, our algorithm reports
lower RMSE values than Particle RAIM. With a fixed number of faults, the
errors generally increase with increasing bias. At a fixed bias value, the
errors decrease with increasing number of faults up to 6 faulty GNSS
measurements since large number of faults are easily excluded by Particle RAIM
producing an improved distribution over the position. The improved
distribution from GNSS further enables the KL divergence metric to exclude
faulty camera measurements and produce a tighter distribution over the
position domain. However, with a higher number of faults, Particle RAIM does
not have enough redundant correct GNSS measurements to exclude the faulty
measurements resulting in higher positioning error. Nevertheless, with added
vision, our algorithm produces better positioning estimates for all test cases
than Particle RAIM.

Figure 6: RMSE from our algorithm and Particle RAIM (baseline) for varying
numbers of faults in GNSS ranging measurements at a fixed added Gaussian bias
value.

Figure 7: RMSE from our algorithm and Particle RAIM (baseline) for various
added Gaussian bias values with fixed number of faulty GNSS measurements.

###  8.2 Integrity Monitoring

We evaluate the integrity risk bounding performance for two Alert Limits, 8 m
and 16 m. For an Alert Limit of 8 m, Table 3 shows that the derived integrity
risk satisfies the performance requirement with very low false alarm and
missed identification rates. While the false alarm rates reported are 0 for
all test cases except two and the missed identification rates are always less
than 0.11. Additionally, the integrity risk bound upper bounds the reference
risk with a failure ratio of less than 0.11 and a bound gap of less than 0.4
for all cases. Figures 8 and 9 further support the observation that the
derived risk bound is able to over bound the reference risk with low failure
rate for the same Alert Limit. The few instances when the derived risk bound
fails to upper bound the reference risk occur due to large sudden jumps in the
reference risk that go undetected considering the fixed size of our motion
samples. However, in general, the integrity risk produced from our algorithm
is able to satisfy the desired performance requirement and successfully
overbound the reference risk for an Alert Limit as small as 8 m. This choice
of Alert Limit is allowed because of the low positioning errors that further
enable non-conservative integrity risk bounds.

Table 3: Integrity Risk for Alert Limit of 8 m  Added Bias Value (meter) | No. of Faults | $P_{FA}$ | $P_{MI}$ | Failure Ratio | Failure Error (meter) | Bound Gap
---|---|---|---|---|---|---
100 | 2 | 0 | 0.03 | 0.07 | 7.5 | 0.26
100 | 4 | 0 | 0.04 | 0.04 | 2.3 | 0.25
100 | 6 | 0 | 0.07 | 0.11 | 2.9 | 0.25
100 | 9 | 0.07 | 0.03 | 0.07 | 4.7 | 0.36
200 | 2 | 0 | 0.07 | 0.07 | 3.5 | 0.20
200 | 4 | 0.11 | 0 | 0.04 | 4.8 | 0.40
200 | 6 | 0 | 0 | 0 | - | 0.38
200 | 9 | 0 | 0.07 | 0.04 | 5.4 | 0.36

Figure 8: Reference risk and integrity risk bound with 8 m Alert Limit for
varying numbers of faults and added bias of 100 m in GNSS measurements. The
derived risk bound over bounds the reference risk with less than 0.11 failure
ratio for all test cases.

Figure 9: Reference risk and integrity risk bound with 8 m Alert Limit for
varying numbers of faults and added bias of 200 m in GNSS measurements. The
derived risk bound over bounds the reference risk with less than 0.07 failure
ratio for all test cases.

For an Alert Limit of 16 m, Table 4 shows that the integrity risk satisfies
the integrity performance requirement with 0 false alarm rates. Furthermore,
the missed identification rates are always 0 except for the test case with 9
faults and 100 m added bias. Specifying a larger Alert Limit lowers the risk
associated with the distribution over position since almost all particles from
the perturbed distributions lie within the Alert Limit. Thus, the integrity
risk with a 16 m Alert Limit is reported to be much smaller compared to the
risk obtained with a 8 m Alert Limit as shown in Figures 8 and 9.
Additionally, the derived risk bound produces even lower failure ratio of less
than 0.07 and a tighter bound gap of less than 0.1. Overall, the derived risk
bound over bounds the reference risk for various bias and fault scenarios in
Figures 10 and 11.

Table 4: Integrity Risk for Alert Limit of 16 m Added Bias Value (meter) | No. of Faults | $P_{FA}$ | $P_{MI}$ | Failure Ratio | Failure Error (meter) | Bound Gap
---|---|---|---|---|---|---
100 | 2 | 0 | 0 | 0 | - | 0.10
100 | 4 | 0 | 0 | 0 | - | 0.08
100 | 6 | 0 | 0 | 0.04 | 5.9 | 0.05
100 | 9 | 0 | 0.04 | 0.07 | 9.7 | 0.08
200 | 2 | 0 | 0 | 0.07 | 5.0 | 0.09
200 | 4 | 0 | 0 | 0.07 | 4.2 | 0.07
200 | 6 | 0 | 0 | 0 | 3.6 | 0.06
200 | 9 | 0 | 0 | 0.04 | 3.8 | 0.01

Figure 10: Reference risk and integrity risk bound with 16 m Alert Limit for
varying numbers of faults and added bias of 100 m in GNSS measurements. The
derived risk bound over bounds the reference risk with less than 0.07 failure
ratio for all test cases.

Figure 11: Reference risk and integrity risk bound with 16 m Alert Limit for
varying numbers of faults and added bias of 200 m GNSS measurements. The
derived risk bound over bounds the reference risk with less than 0.07 failure
ratio for all test cases.

##  9 CONCLUSION

In this paper, we presented a framework for joint state estimation and
integrity monitoring for a GNSS-camera fused system using a particle filtering
approach. To quantify the uncertainty in camera measurements, we derived a
probability distribution directly from camera images leveraging a data-driven
approach along with image registration. Furthermore, we designed a metric
based on KL divergence to probabilistically fuse measurements from GNSS and
camera in a fault-tolerant manner. The metric accounts for vision faults and
mitigates the errors that arise due to cross-contamination of measurements
during sensor fusion. We experimentally validated our framework on real-world
data under NLOS errors, added Gaussian bias noise to GNSS measurements, and
added vision faults. Our algorithm reported lower positioning error compared
to Particle RAIM which uses only GNSS measurements. The integrity risk from
our algorithm satisfied the integrity performance requirement for Alert Limits
of 8 m and 16 m with low false alarm and missed identification rates.
Additionally, the derived integrity risk successfully provided an upper bound
to the reference risk with a low failure rate for both Alert Limits, making
our algorithm suitable for practical applications in urban environments.

##  10 ACKNOWLEDGMENT

We express our gratitude to Akshay Shetty, Tara Mina and other members of the
Navigation of Autonomous Vehicles Lab for their feedback on early drafts of
the paper.

## References

  * [1] S. Gupta and G. X. Gao, “Particle raim for integrity monitoring,” Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation, ION GNSS+ 2019, 2019.
  * [2] H. Zhu, “On information and sufficiency,” 04 1997.
  * [3] N. Zhu, J. Marais, D. Bétaille, and M. Berbineau, “Gnss position integrity in urban environments: A review of literature,” IEEE Transactions on Intelligent Transportation Systems, vol. 19, no. 9, pp. 2762–2778, 2018.
  * [4] Y. C. Lee, “Analysis of range and position comparison methods as a means to provide gps integrity in the user receiver,” Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation), vol. 19, no. 9, pp. 1–4, 1986.
  * [5] S. Bhamidipati and G. Gao, “Slam-based integrity monitoring using gps and fish-eye camera,” Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation, ION GNSS+ 2019, pp. 4116–4129, 10 2019.
  * [6] Z. Wang, Y. Wu, and Q. Niu, “Multi-sensor fusion in automated driving: A survey,” IEEE Access, vol. 8, pp. 2847–2868, 2020.
  * [7] J. Rife, “Collaborative vision-integrated pseudorange error removal: Team-estimated differential gnss corrections with no stationary reference receiver,” IEEE Transactions on Intelligent Transportation Systems, vol. 13, no. 1, pp. 15–24, 2012.
  * [8] He Chengyan, Guo Ji, Lu Xiaochun, and Lu Jun, “Multipath performance analysis of gnss navigation signals,” pp. 379–382, 2014.
  * [9] S. M. Steven Miller, X. Zhang, and A. Spanias, Multipath Effects in GPS Receivers: A Primer.  2015\.
  * [10] K. Ali, X. Chen, F. Dovis, D. De Castro, and A. J. Fernández, “Gnss signal multipath error characterization in urban environments using lidar data aiding,” pp. 1–5, 2012.
  * [11] R. E. Kalman, “A new approach to linear filtering and prediction problems,” Transactions of the ASME–Journal of Basic Engineering, vol. 82, no. Series D, pp. 35–45, 1960.
  * [12] X. Wang, N. Cui, and J. Guo, “Information filtering and its application to relative navigation,” Aircraft Engineering and Aerospace Technology, vol. 81, pp. 439–444, 09 2009.
  * [13] L. Fu, J. Zhang, R. Li, X. Cao, and J. Wang, “Vision-aided raim: A new method for gps integrity monitoring in approach and landing phase,” Sensors (Basel, Switzerland), vol. 15, pp. 22854–73, 09 2015.
  * [14] C. Tanil, S. Khanafseh, M. Joerger, and B. Pervan, “Sequential integrity monitoring for kalman filter innovations-based detectors,” Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation, ION GNSS+ 2018, 10 2018.
  * [15] J. Al Hage and M. E. El Najjar, “Improved outdoor localization based on weighted kullback-leibler divergence for measurements diagnosis,” IEEE Intelligent Transportation Systems Magazine, pp. 1–1, 2018.
  * [16] J. Al Hage, P. Xu, and P. Bonnifait, “Bounding localization errors with student’s distributions for road vehicles,” International Technical Symposium on Navigation and timing, 11 2018.
  * [17] N. A. Tmazirte, M. E. E. Najjar, C. Smaili, and D. Pomorski, “Multi-sensor data fusion based on information theory. application to gnss positionning and integrity monitoring,” 15th International Conference on Information Fusion, pp. 743–749, 2012.
  * [18] Z. Gong, P. Liu, Q. Liu, R. Miao, and R. Ying, “Tightly coupled gnss with stereo camera navigation using graph optimization,” Proceedings of the 32nd International Technical Meeting of the Satellite Division of The Institute of Navigation, ION GNSS+ 2018, pp. 3070–3077, 10 2018.
  * [19] A. Mohanty, S. Gupta, and G. X.Gao, “A particle filtering framework for integrity risk of gnss-camera sensor fusion,” Proceedings of the 33 nd International Technical Meeting of the Satellite Division of The Institute of Navigatin, ION GNSS+ 2020, 2020.
  * [20] P. Reisdorf, T. Pfeifer, J. Breßler, S. Bauer, P. Weissig, S. Lange, G. Wanielik, and P. Protzel, “The problem of comparable gnss results – an approach for a uniform dataset with low-cost and reference data,” in The Fifth International Conference on Advances in Vehicular Systems, Technologies and Applications (M. Ullmann and K. El-Khatib, eds.), vol. 5, p. 8, nov 2016.  ISSN: 2327-2058.
  * [21] F. Gustafsson, F. Gunnarsson, N. Bergman, U. Forssell, J. Jansson, R. Karlsson, and P. . Nordlund, “Particle filters for positioning, navigation, and tracking,” IEEE Transactions on Signal Processing, vol. 50, no. 2, pp. 425–437, 2002.
  * [22] S. B. O. Bousquet and G. Lugosi, “Introduction to statistical learning theory,” Advanced Lectures on Machine Learning, vol. 3176, pp. 169–207, 2004.
  * [23] M. Simandl and J. Dunik, “Design of derivative-free smoothers and predictors,” 14th IFAC Symposium on System Identification, Newcastle, Australia, 03 2006.
  * [24] H. W. Sorenson and D. L. Alspach, “Recursive bayesian estimation using gaussian sums,” Automatica, 1971.
  * [25] J. P. Vila and P. Schniter, “Expectation-maximization gaussian-mixture approximate message passing,” IEEE Transactions on Signal Processing, vol. 61, no. 19, pp. 4658–4672, 2013.
  * [26] C. M. Bishop, Pattern recognition and machine learning.  Information science and statistics, New York, NY: Springer, ”2006”.
  * [27] G. Bradski, “The opencv library,” Dr. Dobb’s Journal of Software Tools, 2000\.
  * [28] E. Rublee, V. Rabaud, K. Konolige, and G. Bradski, “Orb: An efficient alternative to sift or surf,” International Conference on Computer Vision, pp. 2564–2571, 2011.
  * [29] D. Lowe, “Distinctive image features from scale-invariant keypoints,” International Journal of Computer Vision, vol. 60, pp. 91–110, 11 2004.
  * [30] H. Bay, T. Tuytelaars, and L. Van Gool, “Surf: Speeded up robust features,” European Conference on Computer Vision, vol. 3951, pp. 404–417, 07 2006\.
  * [31] P. F. Alcantarilla, J. Nuevo, and A. Bartoli, “Fast explicit diffusion for accelerated features in nonlinear scale spaces,” British Machine Vision Conference, 2013.
  * [32] S. P. Lloyd, “Least squares quantization in pcm,” Information Theory, IEEE Transactions, vol. 28.2, pp. 129–137, 1982.
  * [33] D. Nister and H. Stewenius, “Scalable recognition with a vocabulary tree,” IEEE Computer Society Conference on Computer Vision and Pattern Recognition, vol. 2, pp. 2161–2168, 2006.
  * [34] T. Erven and P. Harremoës, “Rényi divergence and kullback-leibler divergence,” Information Theory, IEEE Transactions on, vol. 60, pp. 3797–3820, 2014.
  * [35] L. G. Valiant, “A theory of the learnable,” Commun. ACM, vol. 27, p. 1134–1142, Nov. 1984.
  * [36] H. Pesonen, “A framework for bayesian receiver autonomous integrity monitoring in urban navigation,” NAVIGATION: Journal of the Institute of Navigation, vol. 58, pp. 229–240, 09 2011.