# “Krein” regularization method

M.V. Takook  _APC, UMR $7164$_
_Université de Paris_
_$75205$ Paris, France_<EMAIL_ADDRESS>
###### Abstract.

A new method of quantum field theory regularization is presented, which is
inspired by the Krein space quantization and quantum metric fluctuations. It
is called the “Krein” regularization method in this paper. We directly recover
physical results previously obtained from QFT standard renormalization
procedure. By applying our approach to the quantum field theory in curved
space-time and quantum linear gravity, we prove that there is no need for the
higher derivative of the metrics tensor for the renormalization of the theory
in the loop correction. The advantages of our method concerning the previous
ones are that the quantum linear gravity is renormalizable.

Proposed PACS numbers: 04.62.+v, 98.80.Cq, 12.10.Dm

##  1\. Introduction

The appearance of singularities in QFT is a manifestation of the existence of
anomalies in this theory. To eliminate these anomalies and to obtain finite
results, the normal ordering procedure and the regularization and
renormalization techniques are applied to QFT. Although these techniques are
efficient to eliminate the singularities and explaining the experimental
results with high accuracy, they do not pertain to the framework of quantum
theory. Moreover, these techniques cannot be applied to quantum field theory
in curved space-time and quantum linear gravity either, regardless of
curvature, i.e. they are not renormalizable.

Negative norm states (negative energy solutions of the field equation) were
first considered by Dirac in $1942$ to deal with these anomalies [1]: “The
appearance of divergent integrals with odd $n$-values in Heisenberg and
Pauli’s form of quantum electrodynamics may be ascribed to the unsymmetrical
treatment of positive- and negative-energy photon states.” Then he presented
an interpretation for the negative and positive probability amplitude. Muto et
Inoue in $1950$ showed that Dirac’s proposal of indefinite metric quantization
has failed to eliminate all divergences in the hole theory. Moreover, there is
difficulty concerning the physical interpretation of the negative probability
[2].

In $1950$, Gupta applied the idea of indefinite metric quantization to the QED
for obtaining a covariant formalism [3]. He quantized the radiation field by
introducing four types of photons: two are transverse, one is longitudinal,
and one is scalar. Although scalar photon states have a negative norm and the
quantization is worked out with indefinite metric quantization, scalar photon
states have positive energy. Then Gupta imposed a condition on the quantum
physical states to eliminate the negative norm states. Hence, the physical
states are the positive norm states whereas the negative norm states are
auxiliary and do not need physical interpretation.

Gupta’s theory was accepted since those negative norm states are just
auxiliary and do not appear in the physical space of states. Although Dirac’s
proposal and interpretation were not accepted, indefinite metric quantization
was viewed as legitimate in many approaches. The perturbation theory for
gravity requires higher derivatives in the free action and their presence in
the Lagrangian leads to ghost fields, states with negative norm [4]. A brief
review of physical problems leading to indefinite metric quantization and non-
hermitian Hamiltonians was presented in [5, 6].

In this regard, i.e. QFT singularity as an anomaly of the theory, it has been
repeatedly speculated that quantum gravitational field may remove the
divergences holding in conventional field theories by providing a natural cut-
off associated with the well known fundamental Planck length [7]. The idea of
using a gravitational field for solving the divergence problem of quantum
field theory was introduced by Deser in $1957$ [8]. Finally, it has been
proven that the singularity of the light cone can only be eliminated by using
the quantum metric fluctuation and some singularities remain. For a review see
[9, 10].

The negative energy solutions of the field equations are discarded for
avoiding negative probability states, but then the symmetrical properties of
the field solutions are broken as was mentioned by Dirac. This fact can be
easily seen in the quantization of the massless minimally coupled scalar field
in de Sitter space-time. In the procedure of its quantization, the elimination
of the negative norm states break the de Sitter invariance [11, 12]. This is
because the positive norm states do not form a covariantly closed complete set
of solutions of the field equations and that the negative norm states for the
massless minimally coupled scalar field cannot be eliminated! By comparison,
one should notice that the positive norm states in Minkowski space-time form a
complete set of solutions and the elimination of the negative norm states does
not break the Poincaré invariance.

For obtaining a covariant and causal quantization of the massless minimally
coupled scalar field in de Sitter space, similar to the Gupta-Bleuler
quantization in QED and the ghost fields in quantum non-abelian gauge theory,
the negative norm states must be taken into account in the theory and should
be considered as auxiliary states. They are not permitted to propagate in the
external legs of the Feynman diagrams and can only propagate in the internal
line. While the Krein space quantization is applied in a rigorous mathematical
way to the free field, the description of an interaction field theory in terms
of the Krein-space quantization approach remains an open mathematical question
[13].

In this paper, regardless of the existence of a rigorous mathematical
formalism for Krein space quantization in the interaction case, we construct a
new method of quantum field regularization by using the Krein space
quantization and the quantum metric fluctuation. We keep Dirac’s and Deser’s
ideas of renormalizability properties of negative norm states and quantum
metric fluctuation respectively and combine them with Gupta’s idea that the
negative norm states are viewed as auxiliary states which are not permitted to
propagate in the external legs of the Feynman diagrams.

In the next section, a brief review of the singularity problems of the Feynman
propagator is given and then we show how the singularities of this function
are eliminated in the Krein space quantization with quantum metric
fluctuations included. Some remarks on Krein space quantization for the
interaction case are given in Section $3$. The new method of regularization is
introduced in Section $4$ and the renormalization conditions for the scalar
field are reviewed. Finally, the effect of this new method of regularization
on the QFT in curved space-time and quantum linear gravity is discussed. In
these cases, there is no need to change the Einstein field equations by the
higher derivatives of the metric for the given renormalization of the theory
in the one-loop approximation.

##  2\. Green function

It is now a well-established fact that the mathematical origin of the
divergences in QFT comes from the singular behavior of the Feynman Green
function. They may be ascribed to the definition of a point (zero) and
infinity on the space-time manifold $\left(\infty\equiv\frac{1}{0}\right)$.
However, these concepts are ill-defined for a quantum theory. The Feynman
propagator for the scalar field is [14]:

$G_{F}^{p}(x,x^{\prime})=-\mathrm{i}<0\mid
T\phi_{p}(x)\phi_{p}(x^{\prime})\mid
0>=\int\frac{\mathrm{d}^{4}k}{(2\pi)^{4}}\frac{e^{-\mathrm{i}k\cdot(x-x^{\prime})}}{k^{2}-m^{2}+\mathrm{i}\epsilon}$
(2.1)
$=-\frac{1}{8\pi}\delta(\sigma)+\frac{m^{2}}{8\pi}\theta(\sigma)\frac{J_{1}\left(\sqrt{2m^{2}\sigma}\right)-\mathrm{i}N_{1}\left(\sqrt{2m^{2}\sigma}\right)}{\sqrt{2m^{2}\sigma}}-\frac{\mathrm{i}m^{2}}{4\pi^{2}}\theta(-\sigma)\frac{K_{1}\left(\sqrt{-2m^{2}\sigma}\right)}{\sqrt{-2m^{2}\sigma}},$

where $2\sigma=\eta_{\mu\nu}(x^{\mu}-x^{\prime\mu})(x^{\nu}-x^{\prime\nu})$ is
the square of geodesic distance between two points $x$ and $x^{\prime}$. The
divergences of this function come out from its singularity behaviour at short
relative distances ($x-x^{\prime}\rightarrow 0$), at large relative distances
($x-x^{\prime}\rightarrow\infty$) and for light cone propagation ($\sigma=0$).
In these limits the divergences of Green functions appear in $\delta(\sigma)$,
$N_{1}$ and $K_{1}$.

It was conjectured that the quantum metric fluctuations can smear out the
singularity of the Green functions on the light cone, i.e. $\delta(\sigma)$
[7, 8, 9, 10]. Due to the quantum metric fluctuation a unique and definite
light cone does not exist and then the Dirac delta singularity disappears. A
flat background space-time together with a linearised perturbation
$h_{\mu\nu}$ propagating upon it constitute the basic modality of quantum
metric fluctuations, i.e.
$g_{\mu\nu}=\eta_{\mu\nu}+h_{\mu\nu},\;\;|h|<|\eta|$. In the unperturbed
space-time, the square of the geodesic separation of points $x$ and
$x^{\prime}$ is defined by
$2\sigma_{0}=\eta_{\mu\nu}(x^{\mu}-x^{\prime\mu})(x^{\nu}-x^{\prime\nu})$. In
a general curved space-time and in the presence of the perturbation
$h_{\mu\nu}$, it becomes
$2\sigma=g_{\mu\nu}(x^{\mu}-x^{\prime\mu})(x^{\nu}-x^{\prime\nu})$, where
$\sigma=\sigma_{0}+\sigma_{1}+O(h^{2})$. $\sigma_{1}$ is the first-order shift
in $\sigma$ (an operator in the linear quantum gravity). The impact of the
metric fluctuation, in this semi-classical approach, on the second and third
parts of the Green function (2.1) is negligible. Its impact on the first term,
which has a delta-function singularity on the light cone, however, could not
be ignored and should be treated as follows. For simplicity, the massless case
is only discussed here.

In the presence of a classical metric perturbation, the Green function for the
massless scalar field has its delta function singularity on the perturbed
light cone, where $\sigma=0$. In general, it may also become nonzero on the
interior of the light cone due to backscattering off of the curvature. We are
primarily interested in the behavior of functions near the perturbed light
cone:

(2.2)
$G_{F}(x-x^{\prime}:m=0)=\frac{-1}{8\pi}\delta(\sigma)=\frac{-1}{16\pi^{2}}\int_{-\infty}^{\infty}\mathrm{d}\alpha
e^{\mathrm{i}\alpha\sigma_{0}}e^{\mathrm{i}\alpha\sigma_{1}}.$

We now replace the classical metric perturbations by a quantum operator in a
vacuum state $|\psi\rangle$. Then $\sigma_{1}$ becomes an operator which is
linear in the gravitational field operator, $h_{\mu\nu}$. $\sigma_{1}$ can be
decomposed into positive and negative frequency parts,
$\sigma_{1}=\sigma^{+}_{1}+\sigma^{-}_{1}$. It can be defined as
$\sigma^{+}_{1}$ and $\sigma^{-}_{1}$ so that $\sigma^{+}_{1}|\psi\rangle=0$
and $\langle\psi\mid\sigma^{-}_{1}=0$. Thus when we average over the metric
fluctuations, the Green function is replaced by its quantum expectation value
[9, 10]:

(2.3) $\langle
G_{F}(x-x^{\prime};m=0)\rangle=\frac{-1}{16\pi^{2}}\int_{-\infty}^{\infty}\mathrm{d}\alpha
e^{\mathrm{i}\alpha\sigma_{0}}e^{-\frac{1}{2}\alpha^{2}\langle\sigma_{1}^{2}\rangle}=\frac{-1}{16\pi^{2}}\sqrt{\frac{\pi}{2\langle\sigma_{1}^{2}\rangle}}\exp\left(-\frac{\sigma_{0}^{2}}{2\langle\sigma_{1}^{2}\rangle}\right).$

This integral converges only if $\langle\sigma_{1}^{2}\rangle>0$. Note that
this averaged Green function is indeed finite at $\sigma_{0}=0$, provided that
$\langle\sigma_{1}^{2}\rangle\neq 0.$ Thus the light cone singularity has been
removed. One can see that in the limit
$\langle\sigma_{1}^{2}\rangle\rightarrow 0$, the Dirac delta function is re-
obtained.

The Krein space quantization was used for a massless minimally coupled scalar
field in de Sitter space-time for preserving the de Sitter invariance and
eliminating the infrared divergence [11, 12]. Moreover, we have shown that the
quantization in Krein space removes also all types of divergences in the Green
function except the light cone singularity. In this procedure, the auxiliary
negative norm states (negative frequency solutions that do not interact with
the physical states or the real physical world) have been utilized. In this
case, the decomposition of the field operator into positive and negative norm
parts reads as:

(2.4) $\phi(x)=\frac{1}{\sqrt{2}}\left[\phi_{p}(x)+\phi_{n}(x)\right]\,,$

where the two parts commute with each other. The “Feynman” propagator or the
time-ordered product propagator in Krein space is [15]:

(2.5) $G_{T}(x,x^{\prime})=-\mathrm{i}<0\mid T\phi(x)\phi(x^{\prime})\mid
0>=\mathrm{Re}\,G_{F}^{p}(x,x^{\prime}).$

It can be rewritten in terms of Dirac delta, Heaviside and Bessel functions
as:

(2.6)
$G_{T}(x,x^{\prime})=\int\frac{\mathrm{d}^{4}k}{(2\pi)^{4}}e^{-\mathrm{i}k\cdot(x-x^{\prime})}\mathcal{P}\frac{1}{k^{2}-m^{2}}=-\frac{1}{8\pi}\delta(\sigma)+\frac{m^{2}}{8\pi}\theta(\sigma)\frac{J_{1}\left(\sqrt{2m^{2}\sigma}\right)}{\sqrt{2m^{2}\sigma}},$

where $\mathcal{P}$ is principal part symbol.

One observes that with Krein space quantization the two-point function has the
light-cone singularity only, and the latter can be removed thanks to the
quantum metric fluctuations. Then by considering the quantum field theory in
Krein space quantization with quantum metric fluctuations included, it is
proved that all singular behaviors of the free scalar Green functions are
completely removed:

(2.7) $\langle
G_{T}(x-x^{\prime})\rangle=-\frac{1}{8\pi}\sqrt{\frac{\pi}{2\langle\sigma_{1}^{2}\rangle}}\exp\left(-\frac{\sigma_{0}^{2}}{2\langle\sigma_{1}^{2}\rangle}\right)+\frac{m^{2}}{8\pi}\theta(\sigma_{0})\frac{J_{1}\left(\sqrt{2m^{2}\sigma_{0}}\right)}{\sqrt{2m^{2}\sigma_{0}}}.$

It is interesting to note that this Green function is causal and real. In the
case of $2\sigma_{0}=0$, due to the quantum metric fluctuation, $h_{\mu\nu}$,
$\langle\sigma_{1}^{2}\rangle\neq 0$, and we obtain:

(2.8) $\langle
G_{T}(0)\rangle=-\frac{1}{8\pi}\sqrt{\frac{\pi}{2\langle\sigma_{1}^{2}\rangle}}+\frac{m^{2}}{8\pi}\frac{1}{2}.$

It should be noted that $\langle\sigma_{1}^{2}\rangle$ is related to the
density of gravitons [9]. In this context the metric fluctuation is defined on
the Minkowski background and the mass $m^{2}$ can be defined concerning this
space-time background.

The Fourier transform of the Green function (2.7) is given by [17]:

(2.9)
$<\tilde{G}_{T}(k)>=\tilde{G}_{1}(k)+\mathcal{P}\frac{m^{2}}{k^{2}(k^{2}-m^{2})}.$

$\tilde{G}_{1}(k)$ is the Fourier transform of the first part of the Green
function. It is a complicated function due to the quantum metric fluctuation
$\langle\sigma_{1}^{2}\rangle$. For the calculation of
$\langle\sigma_{1}^{2}\rangle$ see [10]. In massive field theory for a sake of
simplicity, the propagator $\tilde{G}_{1}(k)$ can be ignored in the one-loop
approximation, and then the second part is only taken into account.

##  3\. Krein space quantization

Let us first remind some features of Krein space quantization in the
interaction case. The Krein space is defined as a direct sum of a Hilbert
space and an anti-Hilbert space. The physical states belong to the Hilbert
space with a positive norm. In general, the effect of the $S$-matrix with
operators acting on the physical subspace of Krein space yields a new state,
which is out of the physical space:

$S|\mbox{Hilbert space}\rangle=|\mbox{Krein space}\rangle\,.$

Then the inner product of this new state with a physical state result in the
projection of this new state onto the Hilbert space as well:

$S_{if}=\langle\mbox{Hilbert space}|S|\mbox{Hilbert
space}\rangle=\langle\mbox{Hilbert space}|\mbox{Krein space}\rangle\,.$

Let us now analyse this $S$-matrix element $S_{if}$. For the sake of
simplicity only the $\lambda\phi^{4}$ theory is discussed. The classical
Lagrangian density and $S$-matrix operator are given respectively by:

(3.1)
$\mathcal{L}_{c}(\phi)=\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi-\frac{1}{2}m^{2}\phi^{2}-\frac{\lambda}{4!}\phi^{4},\;\;\;S=Te^{\mathrm{i}\frac{\lambda}{4!}\int\phi^{4}(x)d^{4}x}.$

Since we are interested to the scattering process the connected diagrams are
first considered. For the two point function at the order $n$ of the
perturbation theory the following integral must be calculated:

$\int d^{4}x_{1}\cdots\int d^{4}x_{n}\langle
q|T\left([\phi(x_{1})]^{4}\cdots[\phi(x_{n})]^{4}\right)|p\rangle^{c}=$ (3.2)
$\int d^{4}x_{1}\cdots\int d^{4}x_{n}\sum_{i_{1}\cdots i_{n}}\langle
q|T\phi(x_{i_{1}})\phi(x_{i_{2}})|p\rangle G_{T}(x_{i_{3}},x_{i_{4}})\cdots
G_{T}(x_{i_{4n-1}},x_{i_{4n}})\,,$

where the sum is over all possible permitted choices for connected diagrams
and $|p\rangle$ and $|q\rangle$ are physical states.

For the calculation of the time-ordered product, the field operators in this
expression must be divided into two types: 1- those which are connected to the
external line (physical states) and 2- those which are connected to the other
field operators. The second type simply produces the time-ordered product
Green function (2.6). For the first type the effect of the negative norm state
cancels out the singularity of the positive norm states, then the theory is
automatically regularized and the normal ordering procedure is not needed. One
can simply prove that the effect of the negative norm states on the first part
of the equation (3.2) in the one-loop approximation, i.e. $\langle
q|T\phi(x)\phi(x)|p\rangle$, cancels out because the unphysical states have
negative norm and are orthogonal to the positive norm states. Hence the usual
result for this part is obtained without using the normal ordering procedure.

It is important to note that the delta function singularity exists in the
second type of products of field operators and that the Krein space
quantization does not remove all of the singularity, as was mentioned by Muto
et Inoue [2]. Therefore the effect of the negative norm states on the internal
line is simply the replacement of the Feynman Green function (2.1) with the
time-ordered product Green function (2.6) and for the external legs, there is
the disappearance of the normal ordering procedure.

Although negative norm states appear in Krein space quantization,
alternatively by imposing two conditions, the negative norm states disappear
from the $S$-matrix elements. The first condition is the “reality condition”
in which the negative norm states do not appear in the external legs of the
Feynman diagram. This condition guarantees that the negative norm states only
appear in the internal legs and the disconnected parts of the Feynman diagram.
The second condition is that the $S$-matrix elements must be renormalized in
the following way:

$S_{if}^{\prime}\equiv\mbox{probability amplitude}=\frac{\langle\mbox{physical
states},in|\mbox{physical states},out\rangle}{\langle 0,in|0,out\rangle}.$

This condition eliminates the negative norm states in the disconnected parts
as well. In the one-loop approximation, one can prove that the negative norm
states disappear and the usual result is obtained [16, 17, 18]. Since the
above proof at all orders of the perturbation theory is not yet completed we
change the perspective and use the Krein space quantization including quantum
metric fluctuation as a new method of quantum field regularisation whereas
this property holds at all orders of perturbation theory.

##  4\. Krein Regularization method

In QFT the S-matrix elements or the probability transition amplitude by using
the LSZ reduction formula, time evolution operator, and Wick’s theorem can be
written in terms of a summation and multiplications of the Feynman Green
functions (2.1). The singularity appears due to the coincident points and also
the multiplication of the Feynman Green functions:

$G_{F}^{p}(x_{i},x_{i})\,,\;\;\left[G_{F}^{p}(x_{i},x_{j})\right]^{n}\,,\;\;\cdots\,.$

Before presenting our regularization method we recall that the regularization
method in QFT is not unique. However, after imposing the renormalization
conditions, a unique result must be obtained for the physical quantities.

By using the above material, we introduce the Krein regularization method. The
procedure can be completed in two simple steps:

  * (a)

Replacing the Feynman Green functions (2.1) with the time-ordered product
propagator (2.7).

  * (b)

Using the same renormalization conditions as the usual method.

From step (a), it is clear that the theory is completely finite and there is
no appearance of any singularity since the two-point function is finite and
free of any divergences. The theory resulting from this replacement becomes
automatically regularized at all orders of the perturbation theory. Therefore
by applying this method to the scalar effective action in curved space-time
one can see that there is no need for higher derivative terms for eliminating
the divergences [19, 20]. Hence following this method makes the quantum linear
gravity renormalizable.

Using step (a) the $S$-matrix element $S_{if}$ is rewritten as:
$\langle\alpha,in|\beta,out\rangle_{reg}=f(\alpha,\beta,\langle\sigma_{1}^{2}\rangle),$
where $f$ is a regular function. The other step (b) guarantees that the
physical result does not change. The effect of the negative norm states and
quantum metric fluctuations for the internal lines is the elimination of the
singularity. The density of gravitons $\langle\sigma_{1}^{2}\rangle$ may be
considered as a regularization parameter. Then, by using the renormalization
conditions, the effect of regularization can be absorbed in the redefinitions
of the physical parameters and we obtain the usual result.

For the massive quantum field theory and in the one-loop approximation one
simply replaces the Feynman Green functions (2.1) with the time-ordered
product propagator (2.6). Then the light cone singularity disappears from the
integral representation, which exists in the probability amplitude.

Finally let us recall the renormalization conditions for the scalar effective
potential or effective action, which is the generating function of the one-
particle irreducible diagrams. The effective Lagrangian or “quantum
Lagrangian” can be established through the loop correction to the classical
Lagrangian:
$\mathcal{L}_{q}=\mathcal{L}_{c}+\hbar\mathcal{L}_{1}+\hbar^{2}\mathcal{L}_{2}+\cdots\,$.
The effect of the negative norm states and the density of gravitons, which
appear in the loop corrections, can be absorbed in the renormalization
procedure. The renormalization conditions are imposed through:

(4.1)
$\lambda_{\mu}=-\left.\frac{\delta^{4}\mathcal{L}_{q}}{\delta\phi^{4}}\right|_{\phi=\mu}\,,\;\;\;m_{\mu}^{2}=-\left.\frac{\delta^{2}\mathcal{L}_{q}}{\delta\phi^{2}}\right|_{\phi=\mu}\,,\;\;\;1=\left.\frac{\delta}{\delta\partial^{\nu}\phi}\frac{\delta\mathcal{L}_{q}}{\delta\partial_{\nu}\phi}\right|_{\phi=\mu}\,.$

The parameter $\mu$ is the energy scale according to which the mass and the
coupling constant are measured.

The effective action in the one-loop approximation in the Krein regularization
approach was previously calculated [16]. $\lambda_{\mu}$ is a function of a
parameter $\mu$ which is the energy scale of the interaction [16]:

$\lambda_{\mu}=\lambda-\frac{\lambda^{2}}{(8\pi)^{2}}\left[6\ln\frac{\mu^{2}}{m^{2}}+19+12\ln
2\right]+O(\lambda^{3}).$

For calculating the running coupling constant a scale of energy must be chosen
as $\phi=\mu=e^{-t}$. Then the running coupling constant is defined as
$\bar{\lambda}(t,\lambda)$. The Beta function can be calculated as well [16]:

(4.2)
$\beta=\frac{d\bar{\lambda}(t,\lambda)}{dt}=\frac{3\lambda^{2}}{16\pi^{2}}\,,$

which is the same as in the usual results in the one-loop approximation.

The QED effective action has been calculated by using the Krein space
quantization and the quantum metric fluctuations at the one-loop approximation
[17]. The radiative corrections of QED have also been calculated using the
Krein regularization method in this approximation [18]. Furthermore, we have
shown that the negative norm states can be used as regularization devices for
certain problems [15, 21, 22, 23, 24].

##  5\. Conclusions

The Krein space quantization for an interaction field has not yet a
satisfactory mathematical construction. In this paper, we have changed the
perspective and the Krein space quantization included quantum metric
fluctuation is considered as a new method of field regularization. The
procedure of this method has been completed in two steps. The first step (a)
not only guarantees that the divergences disappear in the physical quantities
at all orders of perturbation but also guarantees that we do not need any new
terms for the Lagrangian density. This point has great importance for QFT in
curved space-time and quantum linear gravity. The theory is renormalized in
the second step (b) so that the usual results are obtained. The one-loop
effective action for the scalar field in curved space-time has been calculated
with this method [19] and it is proven that one need not change the Einstein
field equations in absorbing the singularity of the scalar effective action,
contrarily to usual previous methods, where the higher derivatives of metric
appear. Our method can be used to calculate physical observables in scenarios
where the effect of quantum linear gravity cannot be ignored. In this model,
the problem of non-renormalizability of linear quantum gravity can be solved
and the theory can be regularized at all order of perturbation theory by
following the first step (a).

Acknowledgements: The author is grateful to Jean-Pierre Gazeau and Eric Huguet
for helpful discussions. The author would like to thank the Collège de France
and l’Université de Paris for their financial support and hospitality.

## References

  * [1] P.A.M. Dirac, Proc. Roy. Soc. A 180, 1 (1942) The physical interpretation of quantum mechanics.
  * [2] T. Muto and K. Inoue, Progress of Theor. Phys. 5, 1033 (1950) Notes on Dirac’s New Quantization Method in the Field Theory.
  * [3] S.N. Gupta, Proc. Phys. Soc. Sect. A 63, 681 (1950) Theory of longitudinal photons in quantum electrodynamics.
  * [4] S.W. Hawking and T. Hertog, Phys. Rev. D 65, 103515 (2002) Living with ghosts.
  * [5] F. Bagarello, J.P. Gazeau, F.H. Szafraniec, M. Znogil, John Wiley and Sons (2015) Non-selfadjoint operators in quantum physics: Mathematical Aspects.
  * [6] A. Ramirez, B. Mielnik, Rev. Mex. Fis. 49S2, 130 (2003) The challenge of non-hermitian structures in physics.
  * [7] B.S. DeWitt, Phys. Rev. Lett. 13, 114 (1964) Gravity: a universal regulator?
  * [8] S. Deser, Rev. Mod. Phys. 29, 417 (1957) General relativity and the divergence problem in quantum field theory.
  * [9] H.L. Ford, Int. J. Theor. Phys. 38, 2941 (1999) Space-time Metric and Light-cone Fluctuations.
  * [10] H. Yu, H.L. Ford, Phys. Rev. D 60, 084023 (1999) Lightcone fluctuations in flat spacetimes with nontrivial topology.
  * [11] S. De Bièvre, J. Renaud, Phys. Rev. D 57, 6230 (1998) The massless quantum field on the 1+1-dimensional de Sitter space.
  * [12] J.P. Gazeau, J. Renaud, M.V. Takook, Class. Quant. Grav. 17, 1415 (2000) Gupta-Bleuler quantization for minimally coupled scalar fields in de Sitter space.
  * [13] T. Garidi, E. Huguet, J. Renaud, J. Phys. A 38, 245 (2005) Krein space quantization in curved and flat space-times.
  * [14] S.A. Fulling, Cambridge Univ. Press (1989) Aspects of Quantum Field Theory in Curved Space-Time.
  * [15] M.V. Takook, Int. J. Mod. Phys. E 11, 509 (2002) Negative norm states in de Sitter space and QFT without renormalization procedure.
  * [16] A. Refaei and M.V. Takook, Mod. Phys. Lett. A 26, 31 (2011) Scalar effective action in Krein space quantization.
  * [17] A. Refaei and M.V. Takook, Phys. Lett. B 704, 326 (2011) QED effective action in Krein space quantization.
  * [18] B. Forghan, M.V. Takook, A. Zarei, Annals of Physics 327, 2388 (2012) Krein regularization of QED.
  * [19] M.V. Takook, Int. J. Mod. Phys. E 14, 219 (2005) A natural renormalization of the one-loop effective action for scalar field in curved space-time.
  * [20] M. Enayati, M.V. Takook, S. Rouhani, Int. J. Theor. Phys. 55, 5055 (2016), Quantum Linear Gravity in de Sitter Universe on Gupta-Bleuler Vacuum State.
  * [21] S. Rouhani and M.V. Takook, EuroPhys. lett. 68, 15 (2004) Tree-level scattering amplitude in de Sitter space.
  * [22] H. Khosravi, M. Naseri, S. Rouhani and M.V. Takook, Phys. Lett. B 640, 48 (2006) Casimir effect in Krein space quantization.
  * [23] H. Pejhan, M.R. Tanhayi, M.V. Takook, Annals of Physics 341, 195 (2014) Casimir effect for a scalar field via Krein quantization.
  * [24] B. Forghan, M. R. Tanhayi, Mod. Phys. Lett. A 30 , 1550126 (2015) Krein regularization of quantum chromodynamics.


$g\left(\ss_n\right)$ commutes with $g\left(\bs{(\overline{1},\overline{m})}\right)$, $g\left(\bs{(\overline{i},\overline{j+1})}\right)$, and $g\left(\bs{(\overline{i+1},\overline{j})}\right)$.

The fact that $g\left(\ss_n\right) = \ss_n$ commutes with $g\left(\bs{(\overline{1},\overline{m})}\right)$ is done in Lemma <ref>(1). We have that
$g\left(\bs{(\overline{i},\overline{j+1})}\right) = \ss_1^{\ss_{m+2}\ss_{m+3}\dotsc \ss_{j+1}\ss_m^{-1}\ss_{m-1}^{-1}\dotsc \ss_{i+1}^{-1}}$, and $g\left(\bs{(\overline{i+1},\overline{j})}\right) = \ss_1^{\ss_{m+2}\ss_{m+3}\dotsc \ss_j\ss_m^{-1}\ss_{m-1}^{-1}\dotsc \ss_{i+2}^{-1}}$ that both obviously commute with $\ss_n$.

(2) Where $i \neq m$ and $j = n-2$, we need to check that
$g\left(\ss_n\right) = \ss_n$ commutes with $g\left(\bs{(\overline{1},\overline{m})}\right)$, $g\left(\bs{(\overline{i+1},\overline{n-2})}\right)$, and check that $g\left(\bs{(\overline{i},\overline{n-1})}\right)g\left(\ss_n\right)g\left(\bs{(\overline{i},\overline{n-1})}\right) = g\left(\ss_n\right) g\left(\bs{(\overline{i},\overline{n-1})}\right) g\left(\ss_n\right)$.

The fact that $\ss_n$ commutes with $g\left(\bs{(\overline{1},\overline{m})}\right)$ and $g\left(\bs{(\overline{i+1},\overline{n-2})}\right)$ is done in Lemma <ref>(1).
For the last relation, we have that

$g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n g\left(\bs{(\overline{i},\overline{n-1})}\right)$     $=$ $\ss_{n-1}^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)\underline{\ss_{n-1} \ss_n \ss_{n-1}^{-1}}g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_{n-1}$

$=$ $\ss_{n-1}^{-1}\underline{g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_n^{-1}}\ss_{n-1}\underline{\ss_ng\left(\bs{(\overline{i},\overline{n-2})}\right)}\ss_{n-1}$

$=$ $\ss_{n-1}^{-1}\ss_n^{-1}\underline{g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_{n-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)}\ss_n\ss_{n-1}$

$=$ $\underline{\ss_{n-1}^{-1}\ss_n^{-1}\ss_{n-1}}g\left(\bs{(\overline{i},\overline{n-2})}\right)\underline{\ss_{n-1}\ss_n\ss_{n-1}}$

by induction hypothesis

$=$ $\ss_n\ss_{n-1}^{-1}\underline{\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)}\ss_n\ss_{n-1}\ss_n$

$=$ $\ss_n\ss_{n-1}^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_n^{-1}\ss_n\ss_{n-1}\ss_n$

$=$ $\ss_ng\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n$.

(3) Where $i \neq m$ and $j = n-1$, we need to check that $\ss_i g\left(\bs{(\overline{i},\overline{n})}\right) \ss_i = g\left(\bs{(\overline{i},\overline{n})}\right) \ss_i g\left(\bs{(\overline{i},\overline{n})}\right)$ and that $g\left(\bs{(\overline{i},\overline{n})}\right)$ commutes with the images by $g$ of elements that correspond to the reflections in the string from item (3) in Lemma <ref>:

$$s_{m+2} \dotsc s_{n-2}s_{n-1}(\overline{i+1},\overline{n-1})s_{i+2} \dotsc s_{m-1}s_m(\overline{1},\overline{m})s_2 \dotsc s_{i-1}.$$

First, we have

$\ss_ig\left(\bs{(\overline{i},\overline{n})}\right) \ss_i $     $=$     $\underline{\ss_i\ss_n^{-1}}g\left(\bs{(\overline{i},\overline{n-1})}\right) \underline{\ss_n\ss_i}$

$=$     $\ss_n^{-1}\underline{\ss_ig\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_i}\ss_n$

$=$     $\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_ig\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n$ by induction hypothesis

$=$     $\underline{\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n}\ss_i\underline{\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n}$

$=$     $g\left(\bs{(\overline{i},\overline{n})}\right)\ss_ig\left(\bs{(\overline{i},\overline{n})}\right)$.

Since $\ss_2$, $\dotsc$, $\ss_{i-2}$, $\ss_{i-1}$, $\ss_{i+2}$, $\dotsc$, $\ss_{m-1}$, $\ss_m$, $\ss_{m+2}$, $\dotsc$, $\ss_{n-3}$, $\ss_{n-2}$ commute with $\ss_n$, then they commute with $g\left(\bs{(\overline{i},\overline{n})}\right) = \ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n$ by applying the induction hypothesis. Because $g\left(\bs{(\overline{1},\overline{m})}\right)$ commutes with $\ss_n$, we also get that $g\left(\bs{(\overline{i},\overline{n})}\right)$ commutes with $g\left(\bs{(\overline{1},\overline{m})}\right)$ by applying the same argument. Now, we prove that $g\left(\bs{(\overline{i},\overline{n})}\right)$ commutes with $\ss_{n-1}$. Actually, we have

$g\left(\bs{(\overline{i},\overline{n})}\right) \ss_{n-1} $     $=$     $\ss_n^{-1}\ss_{n-1}^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)\underline{\ss_{n-1}\ss_n\ss_{n-1}}$

$=$     $\ss_n^{-1}\ss_{n-1}^{-1}\underline{g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_n}\ss_{n-1}\ss_n$

$=$     $\underline{\ss_n^{-1}\ss_{n-1}^{-1}\ss_n}g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_{n-1}\ss_n$

$=$     $\ss_{n-1}\ss_n^{-1}\underline{\ss_{n-1}^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_{n-1}}\ss_n$

$=$     $\ss_{n-1}\underline{\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n}$

$=$     $\ss_{n-1}g\left(\bs{(\overline{i},\overline{n})}\right)$.

Finally, we show that $g\left(\bs{(\overline{i},\overline{n})}\right)$ commutes with $g\left(\bs{(\overline{i+1},\overline{n-1})}\right)$. We have that  $g\left(\bs{(\overline{i},\overline{n})}\right) g\left(\bs{(\overline{i+1},\overline{n-1})}\right)$ is equal to

$\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n\ss_{n-1}^{-1}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\ss_{n-1} = $

$\ss_n^{-1}\ss_{n-1}^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)\underline{\ss_{n-1}\ss_n\ss_{n-1}^{-1}}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\ss_{n-1} = $

$\ss_n^{-1}\ss_{n-1}^{-1}\underline{g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_n^{-1}}\ss_{n-1}\underline{\ss_ng\left(\bs{(\overline{i+1},\overline{n-2})}\right)}\ss_{n-1} = $

$\underline{\ss_n^{-1}\ss_{n-1}^{-1}\ss_n^{-1}}g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_{n-1}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\ss_n\ss_{n-1} = $

$\ss_{n-1}^{-1}\ss_n^{-1}\underline{\ss_{n-1}^{-1}g\left(\bs{(\overline{i},\overline{n-2})}\right)\ss_{n-1}}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\ss_n\ss_{n-1} = $

$\ss_{n-1}^{-1}\ss_n^{-1}\underline{g\left(\bs{(\overline{i},\overline{n-1})}\right)g\left(\bs{(\overline{i+1},\overline{n-2})}\right)}\ss_n\ss_{n-1} = $ by induction hypothesis

$\ss_{n-1}^{-1}\underline{\ss_n^{-1}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)g\left(\bs{(\overline{i},\overline{n-1})}\right)}\ss_n\ss_{n-1} = $

$\ss_{n-1}^{-1}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\underline{\ss_n^{-1}g\left(\bs{(\overline{i},\overline{n-1})}\right)\ss_n}\ss_{n-1} = $

$\ss_{n-1}^{-1}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\underline{g\left(\bs{(\overline{i},\overline{n})}\right)\ss_{n-1}} = $ by the previous case

$\ss_{n-1}^{-1}g\left(\bs{(\overline{i+1},\overline{n-2})}\right)\ss_{n-1}g\left(\bs{(\overline{i},\overline{n})}\right) = $

Suppose that $w_0=w(\overline{m},\overline{j})$, with $m+1\leq j <n$,
so that we are in the situation of Equation <ref>. Let $\refl_1\refl_2 \dotsc \refl_{n-1}$ be the reduced decomposition of $w_0$ described in Lemma <ref>.
Then we can deduce from the relations of the presentation of $\GG$ given in Theorem <ref> that $g\left(\rrefl_{n-1}\right) g\left(\rrefl_{n-2}\right) g\left(\rrefl_{n-1}\right) = g\left(\rrefl_{n-2}\right) g\left(\rrefl_{n-1}\right)  g\left(\rrefl_{n-2}\right)$, and that  $g\left(\rrefl_{n-1}\right)$ commutes with
each of the elements $g\left(\rrefl_{k}\right)$ with $k<n-2$.

We consider the three cases of Lemma  <ref>.

(1) Where  $j<n-2$, we need to check that $\ss_n$ commutes with
$g\left(\bs{(1,j)}\right)$ and $g\left(\bs{(\overline{m},\overline{j+1})}\right)$. This is readily checked since
$$g\left(\bs{(1,j)}\right) = \ss_j^{\ss_{j-1}^{-1}\ss_{j-2}^{-1} \dotsc \ss_2^{-1}} \hbox{ and } g\left(\bs{(\overline{m},\overline{j+1})}\right) = \ss_1^{\ss_{m+2}\ss_{m+3}\dotsc \ss_{j+1}}.$$

(2) Where $j=n-2$, we need to check that $\ss_n$ commutes with
$g\left(\bs{(1,n-2)}\right)$ and $g\left(\bs{(\overline{m},\overline{n-1})}\right)$.
The first check is straightforward.
For the second,
we have that $g\left(\bs{(\overline{m},\overline{n-1})}\right) = \ss_{n-1}^{-1}g\left(\bs{(\overline{m},\overline{n-2})}\right)\ss_{n-1}$
and one shows that
$$g\left(\bs{(\overline{m},\overline{n-1})}\right)\ss_n g\left(\bs{(\overline{m},\overline{n-1})}\right) = \ss_n g\left(\bs{(\overline{m},\overline{n-1})}\right) \ss_n$$

similarly to the case $g\left(\bs{(\overline{i},\overline{n-1})}\right) \ss_n g\left(\bs{(\overline{i},\overline{n-1})}\right) = \ss_n g\left(\bs{(\overline{i},\overline{n-1})}\right) \ss_n$
in the proof of Lemma <ref>(1).

(3) Where $j = n-1$, we have that $g\left(\bs{(\overline{m},\overline{n})}\right) = \ss_1^{\ss_{m+2}\ss_{m+3}\dotsc \ss_n}
= \ss_n^{-1}g\left(\bs{(\overline{m},\overline{n-1})}\right)\ss_n$.
All the commuting relations between $g\left(\bs{(\overline{m},\overline{n})}\right)$ and
$\ss_{m-1}$, $\ss_{m-2}$, $\dotsc$, $\ss_{2}$ are obvious. We are
left to prove that $g\left(\bs{(\overline{m},\overline{n})}\right) \ss_m g\left(\bs{(\overline{m},\overline{n})}\right) = \ss_m g\left(\bs{(\overline{m},\overline{n})}\right) \ss_m$,
and $g\left(\bs{(\overline{m},\overline{n})}\right)$ commutes with $g\left(\bs{(1,n-1)}\right)$, $\ss_{n-1}$, $\ss_{n-2}$, $\dotsc$, $\ss_{m+2}$. This is done similarly to the proof of Lemma <ref>(3).

The next lemma is readily checked.

Suppose that $w_0=w(\overline{i},\overline{n})$, with $1\leq i \leq m$,
so that we are in the situation of Equation <ref> or Equation <ref>. An identical result to the previous lemmas holds in this situation.

This finishes the situation where $w_0$ is of type I. The next two lemmas are for types II and III, respectively.

Suppose that $w_0=w(i,j)$ or $w_0=w(\overline{i},\overline{j})$ with
$1\leq i <j \leq m$,
so that we are in the situation of one of Equations <ref>-<ref>. Let $\refl_1\refl_2 \dotsc \refl_{n-1}$ be the reduced decomposition of $w_0$ as described in Section <ref>.
Then we can deduce from the relations of the presentation of $\GG$ that $g\left(\rrefl_{n-1}\right) g\left(\rrefl_{n-2}\right) g\left(\rrefl_{n-1}\right) =g\left(\rrefl_{n-2}\right) g\left(\rrefl_{n-1}\right)  g\left(\rrefl_{n-2}\right)$, and that $g\left(\rrefl_{n-1}\right)$ commutes with
each of the elements $g\left(\rrefl_{k}\right)$ with $k<n-2$.

In the case of this lemma, we have that $g\left(\rrefl_{n-1}\right) = \ss_n$ and $g\left(\rrefl_{n-2}\right) = \ss_{n-1}$, so the braid relation $\ss_n \ss_{n-1} \ss_n = \ss_{n-1} \ss_n \ss_{n-1}$ is clearly a consequence of the relations of $\GG$. The commuting relations are also obvious since in the situation of Equations <ref>-<ref>, the indices $i,j$ are such that $1 \leq i < j \leq m$, so that they are far away from $n$ (we have $i,j < n-2$). Hence $\ss_n$ obviously commutes with the image by $g$ of the elements $\rrefl_i$ corresponding to the reflections in the reduced decomposition of $w_0$ described in Section <ref>.

Suppose that $w_0=w(i,j)$ or $w_0=w(\overline{i},\overline{j})$ with
$m+1\leq i <j \leq n$,
so that we are in the situation of one of Equations <ref>–<ref>. Let $\refl_1\refl_2 \dotsc \refl_{n-1}$ be the reduced decomposition of $w_0$ as described in Section <ref>.
Then we can deduce from the relations of the presentation of $\GG$ that $g\left(\rrefl_{n-1}\right) g\left(\rrefl_{n-2}\right) g\left(\rrefl_{n-1}\right) =g\left(\rrefl_{n-2}\right) g\left(\rrefl_{n-1}\right)  g\left(\rrefl_{n-2}\right)$, and that $g\left(\rrefl_{n-1}\right)$ commutes with
each of the elements $g\left(\rrefl_{k}\right)$ with $k<n-2$, except for Equation <ref> (with $i=n-1$) where we need to show one additional non-commuting relation.

The argument of the proof is identical to the situation of Equations <ref>–<ref> treated in Lemmas <ref>–<ref>, which is appropriate to leave as an exercise.

Barbara Baumeister, Fakultät für Mathematik, Universität Bielefeld, Postfach 10 01 31, 33501 Bielefeld, Germany,
E-mail<EMAIL_ADDRESS>
Georges Neaime, Fakultät für Mathematik, Universität Bielefeld, Postfach 10 01 31, 33615 Bielefeld, Germany,
E-mail<EMAIL_ADDRESS>
Sarah Rees, School of Mathematics and Statistics, University of Newcastle, Newcastle NE1 7RU, UK,
E-mail<EMAIL_ADDRESS>


# Error State Extended Kalman Filter Multi-Sensor Fusion for Unmanned Aerial
Vehicle Localization in GPS and Magnetometer Denied Indoor Environments

Lovro Marković, Marin Kovač, Robert Milijas, Marko Car, Stjepan Bogdan Authors
are with Faculty of Electrical and Computer Engineering, University of Zagreb,
10000 Zagreb, Croatia (lovro.markovic, marin.kovac, robert.milijas, marko.car,
matko.orsag, stjepan.bogdan) at fer.hr

###### Abstract

This paper addresses the issues of unmanned aerial vehicle (UAV) indoor
navigation, specifically in areas where GPS and magnetometer sensor
measurements are unavailable or unreliable. The proposed solution is to use an
error state extended Kalman filter (ES -EKF) in the context of multi-sensor
fusion. Its implementation is adapted to fuse measurements from multiple
sensor sources and the state model is extended to account for sensor drift and
possible calibration inaccuracies. Experimental validation is performed by
fusing IMU data obtained from the PixHawk 2.1 flight controller with pose
measurements from LiDAR Cartographer SLAM, visual odometry provided by the
Intel T265 camera and position measurements from the Pozyx UWB indoor
positioning system. The estimated odometry from ES-EKF is validated against
ground truth data from the Optitrack motion capture system and its use in a
position control loop to stabilize the UAV is demonstrated.

##  I Introduction

Unmanned aerial vehicles (UAVs) are becoming increasingly popular as their
capabilities such as autonomy, computing power, and payload weight increase.
Thanks to their increased payload weight, UAVs can adopt a diverse sensor
suite including Light Detection and Ranging (LiDAR) sensors, depth and
multispectral cameras, GPS and ultra wideband (UWB) receivers, magnetometer
modules, etc. All of these sensors provide important information for
localization and stabilization of the UAV. Depending on where the UAV is
deployed, reliable and consistent localization is an ongoing problem [1]. In
the context of fusing multiple measurements, different sensor combinations may
prove more useful than others. In this paper, the focus is on indoor
environments where GPS receivers cannot provide an accurate position lock and
magnetometers report inconsistent measurements due to possible electromagnetic
interference in the vicinity.

State estimation techniques related to sensor fusion are often formulated as
filtering and optimization problems. The optimization methods [2],[3], are
more precise, but often require a significant amount of computational
resources. Due to the scope of the state estimation field, this paper focuses
on filtering techniques, mostly based on the extended Kalman filter (EKF) or
its variants, in the context of multi-sensor fusion. In [4], a particle filter
(PF) is used in a land navigation system for sensor fusion and a Monte-Carlo
variant of PF is presented in [5],[6] for UAV localization in GPS denied
environments. A sigma-point Kalman filter (SPKF) is used for integrated
navigation purposes and GPS/IMU fusion in [7] and [8] respectively. The
authors in [9] utilize an unscented Kalman filter (UKF) in a self-calibrating
visual-inertial sensor fusion framework and [10] presents a UKF-based fusion
of UWB and inertial measurements for indoor UAV localization. A novel multi-
state constraint Kalman filter (MSCKF) is presented in [11] used for the
fusion of GPS and visual-inertial odometry (VIO) measurements. Although
suboptimal, regular EKF implementations are still widely used techniques for
state estimation in the context of multi-sensor fusion [12],[13],[14]. The
work presented in [14] is outlined as visual and SLAM sensors are incorporated
into the state estimation for indoor UAV localization.

Figure 1: The figure shows a custom built UAV used for indoor localization
experiments. Reflective markers are attached to the top of the UAV to provide
ground truth measurements from the Optitrack motion capture system. An Ouster
OS1-16 high-resolution LiDAR sensor, when used with Cartographer SLAM,
provides pose sensing capabilities. A Pozyx UWB accurate positioning sensor is
attached to the top of the UAV and provides position measurements. Finally, an
Intel RealSense T265 camera provides visual odometry information.

The error state extended Kalman filter (ES-EKF) relies on linear error state
dynamics for optimal prediction and update of the error state covariance. A
comparative analysis between the classic EKF and the ES-EKF performed in
[15],[16],[17],[18] shows that the ES-EKF is robust to a variety of aircraft
maneuvers as well as imperfect tuning of the sensor noise covariance.
Moreover, it ensures numerical stability and allows the quaternions to be
treated in their minimal representation. The error state is generally small in
magnitude which prevents singularities and gimbal lock from occuring. Finally,
the ES-EKF is more suitable for UAV localization and stabilization due to its
higher computational efficiency. An application of the ES-EKF filter in multi-
sensor state estimation for underwater vehicle localization is presented in
[19]. Other notable work includes [20], in which the authors apply an
iterative extended Kalman filter (I-EKF) based on an error state
representation in a multi-sensor fusion framework and demonstrate the approach
for the case of outdoor UAV localization. It also allows automatic estimation
of the sensor transformations and fusing absolute as well as relative
measurements. The work in [21] also utilizes I-EKF for localization in GNSS-
limited environments.
In this paper, a solution for UAV localization in GPS and magnetometer denied
environments is proposed. Due to its robustness and computational stability,
an ES-EKF is used for reliable state estimation in the context of multi-sensor
fusion. The work presented in [22] gives an overview of quaternion kinematics
and provides an intuitive insight into the technique of error state filtering
and serves as a starting point for the approach proposed in this paper. The
filter implementation is further extended to allow for online drift
estimation. Visual and LiDAR SLAM, commonly used as pose sensing methods in
the context of indoor localization, are known to accumulate drift error over
time, hence the motivation to consider this in state estimation.
The main contributions of this work are as follows. An ES-EKF is adapted to
allow fusion of multiple measurement sources (position, orientation and
velocity) for indoor UAV localization not supported by GPS or magnetometer
measurements. The filter state is also extended to allow position and
orientation drift estimation for each sensor. In addition, an arbiter for
sensor measurements is introduced to eliminate outliers before the fusion
process. Finally, the approach is validated in an indoor laboratory using
motion capture cameras, with the on-board sensor suite containing both visual
and LiDAR SLAM measurements as well as an UWB positioning anchor. The
estimated UAV state is also used for stabilization in a position control loop.
The paper is organized as follows. In Sec. II, an overview of nominal and
error state kinematics is given, along with the equations for the filter
prediction. The approach to fusing measurements from multiple sources using
the ES-EKF formulation is presented in Sec. III. In Sec. IV, an experimental
validation of the described approach is performed in an indoor laboratory
setup with measurements obtained from a variety of sensors suitable for indoor
navigation. Finally, the paper is concluded in Sec. V.

##  II System Kinematics

In this section, the notions of the nominal and the error state are presented
as core concepts of the error state filter formulation. The equations used for
the prediction step in the ES-EKF are shown. Finally, the true state is given
as a combination of the nominal and the error state.

###  II-A Nominal State

The nominal state kinematics are obtained by integrating high-frequency IMU
data $\textbf{u}_{mv}=\left[\textbf{a}_{mv},\bm{\omega}_{mv}\right]^{T}$
containing measured acceleration and angular velocity values respectively. The
nominal state kinematics are presented as follows:

p $\displaystyle=\textbf{p}+\textbf{v}\Delta
t+\frac{1}{2}\left(\text{R}(\textbf{a}_{mv}-\textbf{a}_{b})+\textbf{g}\right)\Delta
t^{2},$ (1) v
$\displaystyle=\textbf{v}+\left(R(\textbf{a}_{mv}-\textbf{a}_{b})+\textbf{g}\right)\Delta
t,$ (2) q
$\displaystyle=\textbf{q}\cdot\textbf{q}\\{(\bm{\omega}_{mv}-\bm{\omega}_{b})\Delta
t\\},$ (3) $\displaystyle\textbf{a}_{b}$ $\displaystyle=\textbf{a}_{b},$ (4)
$\displaystyle\bm{\omega}_{b}$ $\displaystyle=\bm{\omega}_{b},$ (5) g
$\displaystyle=\textbf{g},$ (6) $\displaystyle\textbf{p}_{i}$
$\displaystyle=\textbf{p}_{i},\quad\forall i\in\text{S}_{drift}$ (7)
$\displaystyle\textbf{q}_{i}$ $\displaystyle=\textbf{q}_{i},\quad\forall
i\in\text{S}_{drift},$ (8)

where q{x} is a quaternion representation of the angle-axis vector x,
$\text{S}_{drift}$ is the set of all available sensors, $\textbf{a}_{b}$ and
$\bm{\omega}_{b}$ are acceleration and angular velocity biases respectively.
Estimated drift transformation for each individual sensor is denoted by
position $\textbf{p}_{i}$ and orientation $\textbf{q}_{i}$.

###  II-B Error State

The errors accumulated in the nominal state due to the omitted noise terms and
model inaccuracies are collected in the error state. It has the following
notation:

$\delta\textbf{x}=[\delta\textbf{p},\delta\textbf{v},\delta\bm{\theta},\delta\textbf{a}_{b},\delta\bm{\omega}_{b},\delta\textbf{g},\delta\textbf{p}_{i},\delta\bm{\theta}_{i}]^{T}.$
(9)

The error state angle-axis notation of the UAV orientation and the sensor
drift orientation, $\delta\bm{\theta}$ and $\delta\bm{\theta}_{i}$
respectively, is used as the minimum rotation representation. This does not
expose the system to the risk of singularities or gimbal-lock, due to the
small magnitude of the error state. The same considerations apply to
explaining the accuracy of the linear nature of the system kinematics of the
error state. Small magnitudes of the error state allow the second-order terms
to be omitted without a significant loss of precision. The model of the error
state kinematics with included system noise w is presented as follows:

$\displaystyle\delta x=\text{F}_{x}\delta x+\text{F}_{w}\textbf{w}.$ (10)
$\displaystyle\text{P}=\text{F}_{x}\text{P}{\text{F}_{x}}^{T}+\text{F}_{w}\text{Q}_{w}{\text{F}_{w}}^{T}.$
(11)

The complete derivation of $\text{F}_{x}$ is found in [22]. Since the error
state, in this paper, is extended to include sensor drift estimation
$\text{F}_{x}$ is also extended as follows:

$\displaystyle\text{F}_{x}=\frac{\partial f}{\partial\delta x}=$
$\displaystyle\medmath{\begin{bmatrix}\text{I}_{3}&\text{I}_{3}\Delta
t&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}\\\
0_{3}&\text{I}_{3}&-\text{R}\cdot\text{A}&-\text{R}\Delta
t&0_{3}&\text{I}_{3}\Delta t&0_{3}&0_{3}\\\
0_{3}&0_{3}&\text{R}^{T}\\{\bm{\omega}\\}&0_{3}&-\text{I}_{3}\Delta
t&0_{3}&0_{3}&0_{3}\\\
0_{3}&0_{3}&0_{3}&\text{I}_{3}&0_{3}&0_{3}&0_{3}&0_{3}\\\
0_{3}&0_{3}&0_{3}&0_{3}&\text{I}_{3}&0_{3}&0_{3}&0_{3}\\\
0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&\text{I}_{3}&0_{3}&0_{3}\\\
0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&\text{I}_{3}&0_{3}\\\
0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&\text{I}_{3}\\\ \end{bmatrix}}$ (12)
$\displaystyle\quad\text{A}=\left[(\textbf{a}_{mv}-\textbf{a}_{b})\right]_{\times}\Delta
t,\quad\bm{\omega}=(\bm{\omega}_{mv}-\bm{\omega}_{b})\Delta t,$ (13)
$\displaystyle\text{F}_{w}=\frac{\partial
f}{\partial\textbf{w}}=\begin{bmatrix}0_{3}&0_{3}&0_{3}&0_{3}\\\
\text{I}_{3}&0_{3}&0_{3}&0_{3}\\\ 0_{3}&\text{I}_{3}&0_{3}&0_{3}\\\
0_{3}&0_{3}&\text{I}_{3}&0_{3}\\\ 0_{3}&0_{3}&0_{3}&\text{I}_{3}\\\
0_{3}&0_{3}&0_{3}&0_{3}\\\ 0_{3}&0_{3}&0_{3}&0_{3}\\\
0_{3}&0_{3}&0_{3}&0_{3}\end{bmatrix},$ (14)
$\displaystyle\medmath{\text{Q}_{w}=\begin{bmatrix}\bm{\sigma}_{v}^{2}\Delta
t^{2}\text{I}_{3}&0_{3}&0_{3}&0_{3}\\\ 0_{3}&\bm{\sigma}_{\theta}^{2}\Delta
t^{2}\text{I}_{3}&0_{3}&0_{3}\\\ 0_{3}&0_{3}&\bm{\sigma}_{a}^{2}\Delta
t^{2}\text{I}_{3}&0_{3}\\\ 0_{3}&0_{3}&0_{3}&\bm{\sigma}_{\omega}^{2}\Delta
t^{2}\text{I}_{3}\end{bmatrix}},$ (15)

where R{x} is a matrix rotation representation of the x angle-axis rotation,
$\bm{\sigma}^{2}_{v}$, $\bm{\sigma}^{2}_{\theta}$, $\bm{\sigma}^{2}_{a}$ and
$\bm{\sigma}^{2}_{\omega}$ are velocity, orientation, acceleration bias and
angular velocity bias variances respectively. Skew-symmetric operator is
defined as:

$[\textbf{x}]_{\times}:\mathbf{R}^{3}\longrightarrow\mathfrak{so}(3).$ (16)

##  III Multi-Sensor Fusion

This section shows how the correction step for the ES-EKF implementation is
designed in the context of multi-sensor fusion. Furthermore, the additions
addressing the sensor drift position and orientation are outlined. Baseline
Kalman filter correction step is presented as follows:

K
$\displaystyle=\text{P}\text{H}^{T}(\text{H}\text{P}\text{H}^{T}+\text{V})^{-1},\quad\text{V}=\bm{\sigma}_{s}^{2}\text{I}$
(17) $\displaystyle\widehat{\delta\textbf{x}}$
$\displaystyle=\text{K}(\textbf{y}_{mv}-h(\widehat{\textbf{x}_{t}})),$ (18) P
$\displaystyle=(\text{I}_{3}-\text{K}\text{H})\text{P},$ (19) $\displaystyle
h(\widehat{\textbf{x}_{t}})$
$\displaystyle=\begin{cases}\text{R}\\{\textbf{q}_{i}\\}\widehat{\textbf{p}_{t}}+\textbf{p}_{i},\quad\text{position
sensor},\\\ \textbf{q}_{i}\cdot\widehat{\textbf{q}_{t}},\quad\text{orientation
sensor},\\\ \widehat{\textbf{x}_{t}},\quad\text{otherwise}.\end{cases}$ (20)

True state $\widehat{\textbf{x}_{t}}$ is the current best estimate of the UAV
state. It is expressed as a combination of the nominal state and error state
as follows:

$\displaystyle\widehat{\textbf{p}_{t}}=\textbf{p}+\widehat{\delta\textbf{p}},$
(21)
$\displaystyle\widehat{\textbf{v}_{t}}=\textbf{v}+\widehat{\delta\textbf{v}},$
(22)
$\displaystyle\widehat{\textbf{q}_{t}}=\textbf{q}\cdot\delta\textbf{q}\\{\widehat{\delta\bm{\theta}}\\},$
(23)
$\displaystyle\widehat{\textbf{a}_{b,t}}=\textbf{a}_{b}+\widehat{\delta\textbf{a}_{b}},$
(24)
$\displaystyle\widehat{\bm{\omega}_{b,t}}=\bm{\omega}_{b},+\widehat{\delta\bm{\omega}_{b}}$
(25)
$\displaystyle\widehat{\textbf{g}_{t}}=\textbf{g}+\widehat{\delta\textbf{g}},$
(26)
$\displaystyle\widehat{\textbf{p}_{i,t}}=\textbf{p}_{i}+\widehat{\delta\textbf{p}_{i}},\quad\forall
i\in\text{S}_{drift}$ (27)
$\displaystyle\widehat{\textbf{q}_{i,t}}=\textbf{q}_{i}\cdot\delta\textbf{q}\\{\widehat{\delta\bm{\theta}_{i}}\\},\quad\forall
i\in\text{S}_{drift},$ (28)

The linearized matrix H is defined as follows:

$\text{H}=\frac{\partial h}{\partial\delta\textbf{x}}=\frac{\partial
h}{\partial\textbf{x}_{t}}\frac{\partial\textbf{x}_{t}}{\partial\delta\textbf{x}}.$
(29)

It assumes different forms depending on the type of measurement received:
position, orientation or velocity. For position measurements with drift
estimation the Jacobian H components are derived as follows:

$\displaystyle\frac{\partial
h}{\partial\textbf{x}_{t}}=\begin{bmatrix}\text{R}\\{\textbf{q}_{i}\\}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&\text{I}_{3}&0_{3}\end{bmatrix},$
(30) $\displaystyle\frac{\partial\textbf{x}_{t}}{\partial\delta\textbf{x}}=$
$\displaystyle\medmath{\begin{bmatrix}\text{I}_{3}&0_{3}&0_{3}&0_{3\times
12}&0_{3}\\\ 0_{3}&\text{I}_{3}&0_{3}&0_{3\times 12}&0_{3}\\\ 0_{4\times
3}&0_{4\times
3}&\frac{\partial(\textbf{q}\cdot\delta\textbf{q}\\{\widehat{\delta\bm{\theta}}\\})}{\partial\delta\bm{\theta}}&0_{4\times
12}&0_{4\times 3}\\\ 0_{12\times 3}&0_{12\times 3}&0_{12\times
3}&\text{I}_{12\times 12}&0_{12\times 3}\\\ 0_{4\times 3}&0_{4\times
3}&0_{4\times 3}&0_{4\times
3}&\frac{\partial(\textbf{q}_{i}\cdot\delta\textbf{q}\\{\delta\bm{\theta}_{i}\\})}{\partial\delta\bm{\theta}_{i}}\end{bmatrix}},$
(31) $\displaystyle
where\quad\frac{\partial(\textbf{q}\cdot\delta\textbf{q}\\{\delta\bm{\theta}\\})}{\partial\delta\bm{\theta}}=\frac{1}{2}\begin{bmatrix}-q_{x}&-q_{y}&-q_{z}\\\
q_{w}&-q_{z}&q_{y}\\\ q_{z}&q_{w}&-q_{x}\\\ -q_{y}&q_{x}&q_{w}\end{bmatrix}.$
(32)

Complete proof of Eq. 32 is found in [22]. When fusing other types of
measurements the difference occurs in Eq. 30. For orientation and velocity
measurements the Jacobian component is defined as:

$\displaystyle\medmath{\frac{\partial
h}{\partial\textbf{x}_{t}}=\begin{bmatrix}0_{3}&0_{3}&\text{I}_{3}&0_{3}&0_{3}&0_{3}&0_{3}&\text{R}^{T}\\{\textbf{q}_{mv}\\}\end{bmatrix}},$
(33) $\displaystyle\frac{\partial
h}{\partial\textbf{x}_{t}}=\begin{bmatrix}0_{3}&\text{I}_{3}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}&0_{3}\end{bmatrix},$
(34)

respectively. Finally, the orientation measurement innovation term is
calculated in its angle-axis representation of as follows:

$\textbf{e}_{\theta}=\left[\text{log}\left(R^{T}\\{\textbf{q}_{i}\cdot\widehat{\textbf{q}_{t}}\\}R\\{\textbf{q}_{mv}\\}\right)\right]^{\vee},$
(35)

where log and vee maps are defined as:

log $\displaystyle:\text{SO}(3)\longrightarrow\mathfrak{so}(3),$ (36)
$\displaystyle\left[R\right]^{\vee}$
$\displaystyle:\mathfrak{so}(3)\longrightarrow\mathbb{R}^{3},$ (37)

respectively.

Sensor fault detection is also addressed prior to the observation fusion. The
state estimation methods behind each sensor are inherently complex and error-
prone. Therefore, to ensure that the filter fuses the healthiest measurements,
an arbiter is implemented to eliminate potential outliers. A component-wise
innovation limit test is introduced for each type of measurement. Since a
single sensor can provide reliable position measurements in the horizontal
plane while the vertical axis could be highly noisy, limits of different
magnitudes may be required for each axis. The initial values for the arbiter
limits are chosen as the maximum in-flight deviation from the mean and are
further empirically tuned.

##  IV Experimental Results

In this section, the indoor testing environment is presented along with the
UAV platform and its components. Furthermore, the sensor data and the state
estimation results are compared and analysed to verify the proposed approach.

(a) Position comparison along the X-axis.

(b) Position comparison along the Y-axis.

(c) Position comparison along the Z-axis.

Figure 2: These figures show a comparison of position estimates along all
three axes between Intel Realsense T265 VIO, Pozyx UWB positioning system,
Cartographer SLAM, ES-EKF multi-sensor fusion and Optitrack motion capture
system. The sensor coordinate frames are aligned wrt. each other with the
estimated drift component removed, as explained in Sec. IV-B.

###  IV-A Hardware Setup

A custom-built UAV platform, shown in Fig. 1 assembled by the Kopterworx
company, serves as a base for indoor localization experiments. The PixHawk 2.1
autopilot with ArduCopter firmware resides on a ProfiCNC/HEXKore high-current
power board. An Intel NUC on-board computer runs the entire UAV flight stack
along with the ES-EKF multi-sensor fusion implementation within the ROS
Melodic environment on Ubuntu Linux 18.04. The on-board computer is connected
to the PixHawk autopilot via serial communication. This enables on-board
control of the UAV using ROS, as well as receiving information about the UAV
and the flight controller’s sensors. IMU measurements obtained from PixHawk
through the MAVRos interface are indispensable pieces of information for the
proposed state estimation system. Taking in consideration that PixHawk
provides three IMU sensors, the risk of failure is negligible.
There are three sensors mounted on the UAV that are used for sensor fusion. An
Ouster OS1-16 high-resolution LiDAR sensor is used in combination with
Cartographer SLAM [23] to provide pose measurements. Previous research [24]
has shown that Cartographer SLAM can be used as a pose sensor for UAV
localization and stabilization in outdoor environments. In addition, LiDAR
SLAM is also capable of producing high quality environment maps that are
useful for obstacle navigation, interaction with the environment etc.
especially when combined with a robust state estimation system to stabilize
the UAV. The second sensor is the Pozyx UWB positioning system, the only
sensor used that does not accumulate drift. The third sensor used is an Intel
RealSense T265 camera that provides VIO information. The visual SLAM runs
directly on the camera, freeing up some of the computing power for the on-
board computer.
The indoor laboratory setup, as seen in the video attachment [25], uses the
Optitrack motion capture system to provide ground truth measurements used for
validation purposes only. Although the paper claims to provide robust
localization in magnetometer denied environments, it is difficult to provide
such an elaborate laboratory setup. Therefore, to mimic such conditions, all
internal and external magnetometers on the UAV have been disabled or removed.

###  IV-B ES-EKF Configuration

The parameter set for the ES-EKF is chosen as follows. Acceleration and
angular velocity bias variances are obtained from the IMU specifications data
sheet. The velocity and orientation process variances are tuned to balance the
delay and noise of the filter estimate. The initial values for the sensor
variance parameters are determined by observing signal-to-noise ratio of the
raw sensor measurements. Each variance vector is further tuned based on a
priori sensor measurement knowledge, e.g. Pozyx and Cartographer position
estimates are known to be noisy along the Z-axis.
Furthermore, each sensor’s initial measurement is considered its origin unless
the the sensor does not measure the orientation, in which case a predefined
rotation is used. Since the filter is able to compensate for transformation
inaccuracies online, small errors in the initial sensor transformations are
irrelevant. For accuracy reasons, all position and orientation measurements
have their respective estimated drift removed from the measurement as follows:

$\displaystyle\textbf{p}_{no\,drift}=\text{R}^{-1}\\{\textbf{q}_{i}\\}\left(\textbf{p}_{mv}-\textbf{p}_{i}\right)$
(38)
$\displaystyle\textbf{q}_{no\,drift}=\textbf{q}_{i}^{-1}\cdot\textbf{q}_{mv}$
(39)

Figure 3: This figure shows position error histograms along each of the axes. The histogram values are an estimate of the underlying probability distribution of the error data. The probability distribution estimate is obtained by comparing data against a known probability density function. |  | RMSE |  |
---|---|---|---|---
Estimate | $x[\text{m}]$ | $y[\text{m}]$ | $z[\text{m}]$ | $yaw[^{\circ}]$ |
Cartographer | 0.18 | 0.08 | 1 | 1.15 |
Pozyx | 0.18 | 0.17 | 0.39 | N/A |
Intel T265 | 0.24 | 0.18 | 0.24 | 0.8 |
ES-EKF | 0.1 | 0.1 | 0.31 | 0.49 |
TABLE I: This table shows the root mean square error (RMSE) values for each of
the position axes and yaw angle. Roll and pitch errors are sufficiently small
for each sensor and are therefore omitted. It is important to note that Pozyx
is not an orientation sensor. Figure 4: Map of the indoor laboratory created
by Cartographer during the experimental validation. The room is 10 m by 7.5 m.
Figure 5: This figure shows the same position comparison data as Fig. 2 but on
a smaller time-scale to better showcase the position estimate dynamics. Figure
6: This figure shows a comparison of yaw estimates between Intel Realsense
T265 visual inertial odometry, Cartographer SLAM, ES-EKF multi-sensor fusion
and Optitrack motion capture system. The data is shown on a smaller time-scale
than 2 in order to better showcase the estimated yaw dynamics.

###  IV-C Discussion

Taking into account all the data presented, namely Fig. 4 and Table I, it is
safe to claim that the ES-EKF based multi-sensor fusion is able to achieve
higher overall accuracy than a single sensor. Although the Intel T265 camera
may be slightly more accurate along the z-axis, the Pozyx and Cartographer
measurements along the same axis are far worse than the ES-EKF estimate.
Altitude spikes and heavy drift can be observed in Fig. 2 c) or Fig. 4 \-
Z-axis error plot. Position controlled flight dynamics can be closely examined
in Fig. 5 and Fig. 6. The estimated UAV state indicates no noise or delay wrt.
the ground truth suggesting stable flight conditions (see [25]).

##  V Conclusion

In this paper an ES-EKF multi-sensor fusion method is proposed for indoor UAV
localization in GPS and magnetometer denied environments. In addition to IMU
data, which is indispensable, an Ouster OS1-16 LiDAR is used in combination
with Cartographer SLAM to provide UAV pose measurements, an Intel RealSense
T265 camera provides visual odometry information and a Pozyx UWB positioning
system provides position measurements only. Since SLAM based navigation
methods are known to accumulate drift over time, the ES-EKF model is extended
to account for drift in sensor position and orientation. This method also
helps with sensor calibration inaccuracies. The proposed state estimation
system is validated in an indoor laboratory setup using the Optitrack motion
capture system as the ground truth. To conclude, an ES-EKF used in the context
of a multi-sensor fusion system is overall more accurate than a single-sensor
state estimate. Stable flight dynamics are also demonstrated using estimated
UAV state in a position control loop.

## ACKNOWLEDGMENT

This research was supported by European Commission Horizon 2020 Programme
through project under G. A. number 820434, named ENergy aware BIM Cloud
Platform in a COst-effective Building REnovation Context - ENCORE [26].
Furthermore, this research was a part of the scientific project Autonomous
System for Assessment and Prediction of infrastructure integrity (ASAP) [27]
financed by the European Union through the European Regional Development Fund-
The Competitiveness and Cohesion Operational Programme (KK.01.1.1.04.0041).

## References

  * [1] M. B. Alatise and G. P. Hancke, “A review on challenges of autonomous mobile robot and sensor fusion methods,” IEEE Access, vol. 8, pp. 39830–39846, 2020.
  * [2] X. Ding, F. Han, T. Yang, Y. Wang, and R. Xiong, “Degeneration-aware localization with arbitrary global-local sensor fusion,” Sensors, vol. 21, p. 4042, June 2021.
  * [3] H. Strasdat, J. Montiel, and A. J. Davison, “Visual SLAM: Why filter?,” Image and Vision Computing, vol. 30, pp. 65–77, Feb. 2012.
  * [4] N. Yang, W. F. Tian, Z. H. Jin, and C. B. Zhang, “Particle filter for sensor fusion in a land vehicle navigation system,” Measurement Science and Technology, vol. 16, pp. 677–681, Feb. 2005.
  * [5] F. J. Perez-Grau, R. Ragel, F. Caballero, A. Viguria, and A. Ollero, “An architecture for robust UAV navigation in GPS-denied areas,” Journal of Field Robotics, vol. 35, pp. 121–145, Oct. 2017.
  * [6] P. Carrasco, F. Cuesta, R. Caballero, F. J. Perez-Grau, and A. Viguria, “Multi-sensor fusion for aerial robots in industrial gnss-denied environments,” Applied Sciences, vol. 11, no. 9, 2021.
  * [7] R. van der Merwe, E. Wan, and S. Julier, “Sigma-point kalman filters for nonlinear estimation and sensor-fusion: Applications to integrated navigation,” in AIAA Guidance, Navigation, and Control Conference and Exhibit, American Institute of Aeronautics and Astronautics, June 2004.
  * [8] J. Wendel, J. Metzger, R. Moenikes, A. Maier, and G. F. Trommer, “A performance comparison of tightly coupled GPS/INS navigation systems based on extended and sigma point kalman filters,” Navigation, vol. 53, pp. 21–31, Mar. 2006.
  * [9] J. Kelly and G. S. Sukhatme, “Visual-inertial sensor fusion: Localization, mapping and sensor-to-sensor self-calibration,” The International Journal of Robotics Research, vol. 30, pp. 56–79, Nov. 2010.
  * [10] W. You, F. Li, L. Liao, and M. Huang, “Data fusion of UWB and IMU based on unscented kalman filter for indoor localization of quadrotor UAV,” IEEE Access, vol. 8, pp. 64971–64981, 2020.
  * [11] W. Lee, K. Eckenhoff, P. Geneva, and G. Huang, “Intermittent GPS-aided VIO: Online initialization and calibration,” in 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, May 2020\.
  * [12] G. Nützi, S. Weiss, D. Scaramuzza, and R. Siegwart, “Fusion of IMU and vision for absolute scale estimation in monocular SLAM,” Journal of Intelligent & Robotic Systems, vol. 61, pp. 287–299, Nov. 2010.
  * [13] S. Weiss, M. W. Achtelik, M. Chli, and R. Siegwart, “Versatile distributed pose estimation and sensor self-calibration for an autonomous MAV,” in 2012 IEEE International Conference on Robotics and Automation, IEEE, May 2012.
  * [14] E. López, R. Barea, A. Gómez, Á. Saltos, L. M. Bergasa, E. J. Molinos, and A. Nemra, “Indoor SLAM for micro aerial vehicles using visual and laser sensor fusion,” in Advances in Intelligent Systems and Computing, pp. 531–542, Springer International Publishing, Dec. 2015.
  * [15] W. Liu, D. Song, Z. Wang, and K. Fang, “Comparative analysis between error-state and full-state error estimation for KF-based IMU/GNSS integration against IMU faults,” Sensors, vol. 19, p. 4912, Nov. 2019\.
  * [16] V. Madyastha, V. Ravindra, S. Mallikarjunan, and A. Goyal, “Extended kalman filter vs. error state kalman filter for aircraft attitude estimation,” in AIAA Guidance, Navigation, and Control Conference, American Institute of Aeronautics and Astronautics, June 2011.
  * [17] Guo, Zhou, Zhu, and Bai, “Low-cost sensors state estimation algorithm for a small hand-launched solar-powered UAV,” Sensors, vol. 19, p. 4627, Oct. 2019.
  * [18] W. Youn and S. Andrew Gadsden, “Combined quaternion-based error state kalman filtering and smooth variable structure filtering for robust attitude estimation,” IEEE Access, vol. 7, pp. 148989–149004, 2019.
  * [19] N. Shaukat, A. Ali, M. J. Iqbal, M. Moinuddin, and P. Otero, “Multi-sensor fusion for underwater vehicle localization by augmentation of RBF neural network and error-state kalman filter,” Sensors, vol. 21, p. 1149, Feb. 2021.
  * [20] S. Lynen, M. W. Achtelik, S. Weiss, M. Chli, and R. Siegwart, “A robust and modular multi-sensor fusion approach applied to MAV navigation,” in 2013 IEEE/RSJ International Conference on Intelligent Robots and Systems, IEEE, Nov. 2013.
  * [21] T. Bulunseechart and P. Smithmaitrie, “A method for uav multi-sensor fusion 3d-localization under degraded or denied gps situation,” Journal of Unmanned Vehicle Systems, vol. 6, no. 3, pp. 155–176, 2018.
  * [22] J. Solà, “Quaternion kinematics for the error-state kalman filter,” 2017.
  * [23] W. Hess, D. Kohler, H. Rapp, and D. Andor, “Real-Time Loop Closure in 2D LIDAR SLAM,” in 2016 IEEE International Conference on Robotics and Automation (ICRA), pp. 1271–1278, 2016.
  * [24] R. Milijas, L. Markovic, A. Ivanovic, F. Petric, and S. Bogdan, “A comparison of lidar-based slam systems for control of unmanned aerial vehicles,” in 2021 International Conference on Unmanned Aircraft Systems (ICUAS), pp. 1148–1154, 2021.
  * [25] L. Markovic and M. Kovac and R. Milijas and M. Car and S.Bogdan, “Error state extended kalman filter multi-sensor fusion forunmanned aerial vehicle localization in gps and magnetometerdenied indoor environments.” https://www.youtube.com/playlist?list=PLC0C6uwoEQ8aA2cqpzYoNYVxpQuZh4YaI.
  * [26] “Encore project.” http://encorebim.eu/.  Accessed: 2019-09-10.
  * [27] “Asap project.” https://asap-project.com/en/.  Accessed: 2021-02-27.