###### Abstract

This work presents a number of techniques to improve the ability to create
magnetic field maps on a UAV which can be used to quickly and reliably gather
magnetic field observations at multiple altitudes in a workspace.
Unfortunately, the electronics on the UAV can introduce their own magnetic
fields, distorting the resultant magnetic field map. We show methods of
reducing and working with UAV-induced noise to better enable magnetic fields
as a sensing modality for indoor navigation. First, some gains in our flight
controller create high-frequency motor commands that introduce large noise in
the measured magnetic field. Next, we implement a common noise reduction
method of distancing the magnetometer from other components on our UAV.
Finally, we introduce what we call a compromise GPR (Gaussian process
regression) map that can be trained on multiple flight tests to learn any
flight-by-flight variations between UAV observation tests. We investigate the
spatial density of observations used to train a GPR map then use the
compromise map to define a consistency test that can indicate whether or not
the magnetometer data and corresponding GPR map are appropriate to use for
state estimation. The interventions we introduce in this work facilitate
indoor position localization of a UAV whose estimates we found to be quite
sensitive to noise generated by the UAV.

###### keywords:

Gaussian process; magnetic field mapping; unmanned aerial vehicle (UAV);
fingerprinting; indoor environments

1 1 0 https://doi.org/ Fast and Noise-Resilient Magnetic Field Mapping on a
Low-Cost UAV Using Gaussian Process Regression Fast and Noise-Resilient
Magnetic Field Mapping on a Low-Cost UAV Using Gaussian Process Regression
Prince E. Kuevor 1,*, Maani Ghaffari2, Ella M. Atkins 3, and James W. Cutler 4
Prince E. Kuevor, James W. Cutler and Ella M. Atkins Kuevor, P.E.; Ghaffari,
M.; Atkins, E.M.; Cutler, J.W. Correspondence<EMAIL_ADDRESS>(P.E.K.)

##  1 Introduction

Earth’s magnetic field provides a reference measurement anywhere on the planet
to assist with navigation. Magnetometers, or 3D compasses, onboard many modern
robotics platforms can measure the local magnetic field to assist with
navigation. To do this, we often rely on simplifying assumptions (i.e., the
magnetic field points northward and is constant in a target region) or on maps
of Earth’s magnetic field like the World Magnetic Model (WMM). However, when
testing in indoor environments, these assumptions and worldwide maps are
inaccurate due to the contribution of the magnetic field from metallic objects
inside buildings. Because of this, measurements from magnetometers are often
ignored for indoor navigation.

There are a number of existing works that leverage indoor magnetic field
measurements to estimate position, or attitude Kuevor et al. (2021); Kok and
Solin (2018); Vallivaara et al. (2011, 2010); Akai and Ozaki (2015, 2017);
Suksakulchai et al. (2000); Li et al. (2012); Haverinen and Kemppainen (2009);
Wu et al. (2019). However, a common theme with the position localization and
SLAM (simultaneous localization and mapping) papers Robertson et al. (2013);
Vallivaara et al. (2010, 2011); Wang et al. (2017) is that they do not
investigate the accuracy their magnetic map. Thus it is not clear if their
maps can be used as a valuable prior for future experiments in the mapped
region. Similarly, for indoor attitude estimation, many previous works see
Earth’s magnetic field (mapped by WMM) as the signal and any magnetic fields
from the building itself to be a source of noise that must be rejected Afzal
et al. (2011); Yin and Zhang (2018); Karimi et al. (2020). Though this
paradigm works for attitude estimates, it strips away valuable information
that is essential to using magnetic fields for position localization since
Earth’s magnetic field is fairly constant over many kilometers.

Figure 1: Methods to reduce, incorporate, and identify vehicle-induced
magnetic field noise. Our suggested interventions can be progressively added
until the user’s GPR map is sufficiently accurate for the application.

This paper, outlined in Figure 1, aims to augment indoor sensor suites with
magnetic field data by demonstrating how to create and validate magnetic field
maps of indoor spaces. The key idea is to enable the creation of local indoor
magnetic field maps that can then be leveraged to improve position and
attitude estimates. The value of indoor magnetic field maps is demonstrated in
works that use them to estimate the position of pedestrians and robotic
vehicles inside buildings and to estimate the attitude of a drone indoors. As
shown by ref. Akai and Ozaki (2017), gathering observations for a magnetic
field map by hand is time-consuming and does not scale to mapping increasingly
large spaces. As such, we propose using an unmanned aerial vehicle (UAV) to
gather observations in the target volume of the workspace, using Gaussian
Process Regression (GPR) to interpolate between these observations, then
finally leveraging the map in future experiments within the mapped volume.

The contributions of this work are as follows.

  1. 1.

The magnetic noise induced by the motors and ESCs on a UAV can be reduced by
removing high-frequency commands from the flight controller to the motors.

  2. 2.

Our UAV produces magnetic biases that occasionally change in magnitude and
direction. We were unable to identify the underlying cause but found that
distancing the magnetometer from the electronics, a commonly used intervention
to reduce the measured magnetic field noise, decreased the measured variation
of these biases.

  3. 3.

Our “compromise map” can be trained on large datasets to prevent overfitting
to a single flight test and learn the flight-by-flight variations of the drone
without incurring high computation costs when predicting the magnetic field.

  4. 4.

We find that our compromise map has similar accuracy if the location of
training observations is no further than 0.55m apart. This agrees with a
similar, qualitative study by Akai and Ozaki in ref. Akai and Ozaki (2017).

  5. 5.

We introduce a consistency test to indicate if a given GPR map has good
agreement with new measurements from a subsequent flight test. This
consistency test was important for contextualizing our position localization
accuracy in another project.

  6. 6.

In a comparison of two methods of mapping the norm (magnitude) of the magnetic
field, we find (empirically) that creating a specialized map on the norm of
the field is equivalent to taking the norm of estimates from a vector-valued
magnetic field map.

At a high level, this paper rings similar to our previous work in Kuevor et
al. (2021) where we used a drone to create a magnetic field map, interpolate
with GPR, then finally demonstrate the value of the map in improving attitude
estimates of the drone. This work focuses on improving the accuracy of
magnetic field maps in the face of vehicle-induced noise and comparing
different methods of creating magnetic field maps.

The remainder of this paper is structured as follows. Section 2 presents
related works, and Section 3 introduces our mathematical notation for Gaussian
process regression and our compromise map. Section 4 introduces our UAV and
testing methodology for this work. Finally, Section 5 presents our results,
while conclusions and future work appear in Section 6.

##  2 Related Works

A magnetic field map must have an input dimension $p$ and output dimension
$m$, which we will denote as a $p\rightarrow m$ map. The choice of input
dimension $p$ often depends on the agent making the map. For example, $p=2$
for wheeled robots Li et al. (2012); Vallivaara et al. (2011, 2010); Wu et al.
(2019); Akai and Ozaki (2015); Hanley et al. (2021) and most pedestrian
localization while $p=3$ for UAVs Brzozowski et al. (2016); Brzozowski and
Kazmierczak (2017), multi-floor pedestrian localization Kok and Solin (2018),
or specially-outfitted ground robots Akai and Ozaki (2017); Hanley et al.
(2021). The output $m$ is set by how the map will be utilized. Works like
Suksakulchai et al. (2000); Haverinen and Kemppainen (2009) use $m=1$ for the
orientation or magnitude of the magnetic field vector, while other works like
Kuevor et al. (2021); Vallivaara et al. (2010, 2011); Wu et al. (2019); Li et
al. (2012); Akai and Ozaki (2015, 2017); Kok and Solin (2018) use all three
vector components of the magnetic field with $m=3$.

A number of works using indoor magnetic field maps are focused on pedestrian
localization. Such works often utilize foot-mounted, or calf-mounted sensor
suites Frassl et al. (2013); Robertson et al. (2013); Liu et al. (2021) that
can leverage zero-velocity updates to help with localization. Alternatively,
many such works utilize the ubiquity of smartphones Gozick et al. (2011) and
even consider constructing maps by crowdsourcing magnetic field measurements
from many users Ayanoglu et al. (2018). To our knowledge, all pedestrian-
localization-based magnetic field maps are $2\rightarrow m$ (or
“$2.5\rightarrow m$” maps Liu et al. (2021)), which assume a constant (or
near-constant) altitude.

Ground vehicles are another common platform for indoor magnetic field mapping
Li et al. (2012); Vallivaara et al. (2011, 2010); Wu et al. (2019); Akai and
Ozaki (2015); Hanley et al. (2021); Almeida et al. (2021); Frassl et al.
(2013). These platforms benefit from wheel encoder odometry to assist their
state estimates. Most of such works also create $2\rightarrow m$ maps since
the platform is confined to the ground. However, some researchers outfit their
wheeled robots with a magnetometer on a vertical actuator Akai and Ozaki
(2017) or vertically-spaced magnetometers Hanley et al. (2021) that let them
create $3\rightarrow m$ maps. The latter is a recent work by Hanley et al.
showing that indoor magnetic fields can be quite sensitive to changes in
altitude Hanley et al. (2021) and emphasizing some issues with $2\rightarrow
m$ planar magnetic field maps.

UAVs are used more regularly for outdoor magnetic field surveys than for
indoor mapping. A recent review paper by Zhang et al. Zheng et al. (2021)
presents works on outdoor magnetic field surveillance and also addresses
methods of characterizing and suppressing UAV-induced magnetic noise. From the
papers in their review, it is clear that the size and cost of the outdoor
survey vehicles and sensors tend to be much larger. For example, Versteeg and
McKay (2007) mounted a Geometrics G823A cesium magnetometer on a gas-powered
helicopter and measured 800nT of magnetic variation caused by components on
the vehicle when it is not powered. This reduced to 80nT and 40nT when the
magnetometer was attached to a boom of length 0.5m and 1.2m, respectively.
Additionally, they consider vehicles that cost $2K - $45K USD and weigh 7lbs -
51lbs. The larger vehicle size of outdoor magnetic surveys allows for
interventions like suspending a magnetometer at the end of a 4.5m-long cable
Koyama et al. (2013). Thus, the vehicle size and project budget of outdoor
surveying techniques allow for solutions to vehicle-induced magnetic noise
(and magnetometers with better sensing capabilities) that are not applicable
to (or not used on) indoor platforms. Zheng et al. cover other UAV noise
characterization and mitigation efforts for outdoor magnetic surveys in
Section 3 of ref. Zheng et al. (2021).

By contrast, some works use UAVs with magnetometers indoors, but do not create
a map of the magnetic field for their implementation. Brzozowski et al.
describe methods to support the use of indoor magnetic fields on UAVs in refs.
Brzozowski et al. (2016); Brzozowski and Kazmierczak (2017), but mostly
present methods of gathering and visualizing magnetic field observations. They
stop short of interpolating their magnetic field observations to create a
queryable map or perform any state estimation. Furthermore, they do not
actually use a flight vehicle (real or simulated) to gather or validate their
magnetic field measurements. The authors in ref. Vasiljevic et al. (2022)
perform 3D position localization of a UAV near a two-wire power line setup in
their lab space, but do not attempt to map the magnetic field near the wires
in the process. Li et al. estimate the 6DOF pose of a UAV indoors using
magnetometers (among other sensors), but only use the magnetometer for heading
estimates and achieve this without a magnetic field map Li et al. (2019).
Zahran et al. leverage hall effect sensors (a type of magnetometer) in a
clever way to estimate the velocity of a quadrotor and improve the dead
reckoning of such flight vehicles (again, without magnetic field maps) Zahran
et al. (2019).

There are works that use UAVs to map signals like ultra-wideband (UWB) radio
Tiemann et al. (2018) and visual light communication (VLC) Niu et al. (2021).
Some works often use the term “fingerprinting” to describe the association of
unique features in a reference signal to specific locations in a workspace. In
this work, we describe the same process as mapping.

To our knowledge, however, there are only two works that use a UAV (drone,
multirotor, flight vehicle, etc…) to make or leverage indoor magnetic field
maps. One is our previous work Kuevor et al. (2021) which uses Gaussian
process regression (GPR) to make $3\rightarrow 3$ magnetic field maps of the
DC magnetic field in an indoor workspace of dimensions
[4$\times$3$\times$2.25]m. The other work, by Lipovsky et al. Lipovský et al.
(2021), uses observations from a UAV and cubic spline interpolation to create
$3\rightarrow 3$ maps of the DC, 50Hz, and 150Hz magnetic fields in a
[5$\times$5$\times$2.5]m room. Ref. Lipovský et al. (2021) mounts their
magnetometer 40cm away from the drone based on lessons learned from their
previous investigation on drone-induced magnetic noise Blažek et al. (2018).

In this work, we use a UAV to create $3\rightarrow 3$ indoor magnetic field
maps. A flight vehicle allows for quick and repeatable mapping and validation
of an indoor volume, but can also impede accurate measurements of the ambient
field with magnetic noise generated by the motors and ESCs (electronic speed
controllers). The focus of this work is on identifying the type of magnetic
noise our UAV creates, reducing the amount of magnetic noise, and finally
incorporating any remaining measurable variations into the magnetic field map.

Additionally, we study how the spatial density of magnetic observations affect
the accuracy of a magnetic field map (similar to a study done in Akai and
Ozaki (2017)) and empirically demonstrate the equivalence of two maps that
estimate the norm of the ambient magnetic field. Finally, we introduce a
consistency check to reason about when a GPR-based magnetic field map is
appropriate to use for state estimation.

The interventions and studies done in this paper are lessons we learned while
investigating the 3D position localization of a UAV using a low-cost IMU and
magnetometer. Though a presentation of our position localization methods and
accuracy is outside the scope of this paper, it is important to emphasize that
some of our interventions may not be needed for other applications of indoor
magnetic field maps (e.g., they were not applied for our study on attitude
estimation using indoor magnetic field maps Kuevor et al. (2021)).

##  3 Mathematical Preliminaries and Methodology

This section introduces our use of Gaussian Process Regression (GPR) to
interpolate a set of magnetic field observations in our workspace. Special
notation is used to distinguish a set of $n_{2}$ observations used to train
hyperparameters and a separate set of $n_{1}$ observations used to perform
inference. Additionally, we introduce performance metrics used later in our
analysis of our GPR-based magnetic field maps.

We first introduce some notation. A single measurement of the magnetic field
at an unspecified location is $\tilde{\bm{y}}\in\mathbb{R}^{3}$ with the $x$,
$y$, and $z$ components of this measurement denoted as $\tilde{y}_{x}$,
$\tilde{y}_{y}$, $\tilde{y}_{z}\in\mathbb{R}$ respectively. In general, a
subscript of x, y, or z denotes that respective component of the magnetic
field while an overhead tilde $\tilde{}$ denotes a measured value.
$\tilde{\bm{Y}}^{n}\in\mathbb{R}^{3\times n}$ is a set of $n$ magnetic field
measurements while $\tilde{\bm{Y}}^{n}_{z}\in\mathbb{R}^{n}$ denotes the $z$
component of each magnetic field measurement in the set. Similarly,
$\bm{X}^{n}\in\mathbb{R}^{3\times n}$ is a collection of $n$ spatial locations
in our workspace. The predicted or estimated magnetic field $\hat{\bm{m}}$ at
some location $\bm{r}\in\mathbb{R}^{3}$ is denoted as
$\hat{\bm{m}}(\bm{r})\in\mathbb{R}^{3}$ and the $x$ component of this
prediction is $\hat{m}_{x}(\bm{r})\in\mathbb{R}$. Thus, a collection of $n$
magnetic field estimates will be
$\hat{\bm{M}}^{n}(\bm{X}^{n})\in\mathbb{R}^{3\times n}$. Similarly,
$\hat{\bm{M}}^{n}$ is a set of predictions of the magnetic field where the
location of these predictions is arbitrary and
$\hat{\bm{M}}^{n}_{y}\in\mathbb{R}^{n}$ gives just the predictions from
$\mathcal{GP}_{y}$. Generally, an overhead hat $\hat{}$ denotes a
prediction/estimate while a superscript integer $n$ denotes the number of
measurements/predictions/locations of the respective matrix.

###  3.1 Gaussian Process Regression

Gaussian Process Regression (GPR) is a machine learning tool that can be used
to estimate a signal given a set of noisy measurements. Here, we use GPR as
the backbone of our magnetic field map leveraging methods from Rasmussen et
al. in Rasmussen and Williams (2006). The goal is to have a magnetic field map
of the flight workspace which will provide an estimated magnetic field vector
at any location in the working volume.

This work uses GPR to create three separate $3\rightarrow 1$ maps each
responsible for estimating the $x$, $y$, and $z$ components of the magnetic
field vector in our workspace. To create the GPR-based map, we first gather
$n$ observations of the magnetic field throughout the workspace. Observation
sets give the three components of the measured magnetic field
$\tilde{\bm{Y}}^{n}_{x}$, $\tilde{\bm{Y}}^{n}_{y}$,
$\tilde{\bm{Y}}^{n}_{z}\in\mathbb{R}^{n}$ at each 3D position in our design
matrix $\bm{X}^{n}\in\mathbb{R}^{3\times n}$. Together, these quantities
define the training sets $D_{x}=(\bm{X},\tilde{\bm{Y}}_{x})$,
$D_{y}=(\bm{X},\tilde{\bm{Y}}_{y})$, and $D_{z}=(\bm{X},\tilde{\bm{Y}}_{z})$
for the $x$, $y$, and $z$ GPRs respectively.

In Ref. Rasmussen and Williams (2006), Rasmussen and Williams define a
Gaussian process as a distribution over functions written as

$f(\bm{r})=\mathcal{N}\left(\bar{f},\mathbb{V}(f)\right)\sim\mathcal{GP}\left(\bm{0},k(\bm{r},\bm{r}^{\prime})\right)$
(1)

where $\bm{0}$ is the zero-mean function, $k$ is the covariance function (or
kernel) and $\bm{r}$, $\bm{r}^{\prime}\in\mathbb{R}^{3}$ are 3D positions.
This work uses the squared exponential covariance function

$k(\bm{r},\bm{r}^{\prime})=\sigma_{f}^{2}\exp\left(-\frac{1}{2l^{2}}(\bm{r}-\bm{r}^{\prime})^{\top}(\bm{r}-\bm{r}^{\prime})\right)+\sigma_{n}^{2}$
(2)

with $\sigma_{f}$ as the signal variance, $l$ as the length scale, and
$\sigma_{n}$ as the sensor noise variance which define our set of
hyperparameters $\bm{\Theta}=\\{\sigma_{f},l,\sigma_{n}\\}$ for this kernel.
In essence, the squared exponential kernel is a measure of similarity between
two 3D locations $\bm{r}$ and $\bm{r}^{\prime}$ (scaled by the
hyperparameters) and will serve as a weighting term in the inference of the
magnetic field at unobserved locations.

Recall that we perform regression using three Gaussian processes
$\mathcal{GP}_{x}$, $\mathcal{GP}_{y}$, and $\mathcal{GP}_{z}$ to represent
the full magnetic field anywhere in the workspace. As such, we need three sets
of hyperparameters ($\bm{\Theta}_{x}$, $\bm{\Theta}_{y}$, $\bm{\Theta}_{z}$)
each computed using the Gaussian Processes for Machine Learning gpml-matlab
toolbox by minimizing the negative log marginal likelihood of the respective
training sets $D_{x}$, $D_{y}$, and $D_{z}$. We define hyperparameters
optimized over the observations in a training set $D_{x}$ as
$\bm{\Theta}^{*}_{x}(D_{x})$. (gpml-matlab was created by Carl Edward
Rasmussen and Hannes Nickisch
http://www.gaussianprocess.org/gpml/code/matlab/doc/ accessed on September
2021).

Now, to estimate the magnetic field at some location $\bm{r}^{*}$, we need a
training set $D$ and a set of hyperparameters $\bm{\Theta}$ corresponding to a
selected kernel $k$. A squared exponential kernel $k$ using optimal
hyperparameters $\bm{\Theta}^{*}(D)$ is denoted as
$k_{\bm{\Theta}^{*}(D)}(\bm{r},\bm{r}^{*})$ where the subscript on the kernel
$k$ specifies which set of hyperparameters are used and $\bm{r}\in\bm{X}$ is a
location in the training set.

This is the extent of our notation, but an example is helpful. Say we want to
estimate the $z$ component of the magnetic field $\hat{m}_{z}(\bm{r}^{*})$ at
some location $\bm{r}^{*}$. For this, $\mathcal{GP}_{z}$ requires a kernel $k$
to compare the target location $\bm{r}^{*}$ against locations $\bm{X}^{n_{1}}$
in its training set
$D_{z}^{n_{1}}=(\bm{X}^{n_{1}},\tilde{\bm{Y}}_{z}^{n_{1}})$. Additionally, its
hyperparameters $\bm{\Theta}^{*}(D_{z}^{n_{2}})$ are optimized over a separate
set of observations
$D_{z}^{n_{2}}=(\bm{X}^{n_{2}},\tilde{\bm{Y}}_{z}^{n_{2}})$. Mathematically,
we express this as

$\begin{split}\hat{m}_{z}(\bm{r}^{*})&=GPR_{z}(\bm{r}^{*})\\\
&=\mathbb{E}\left[\mathcal{GP}_{z}\left(\bm{0},\;K_{\bm{\Theta}_{z}^{*}(D_{z}^{n_{2}})}(\bm{X}^{n_{1}},\bm{r}^{*})\right)\right]\end{split}$
(3)

where a matrix $K(\bm{X},\bm{r}^{*})$ has scalar elements
$k\left(\bm{X}^{\\{i\\}},\bm{r}^{*}\right)$ with each
$\bm{X}^{\\{i\\}}\in\mathbb{R}^{3}$ a column of $\bm{X}\in\mathbb{R}^{3\times
n}$ while $\hat{m}_{z}(\bm{r}^{*})\in\mathbb{R}$. Finally,
$GPR_{z}(\bm{r}^{*})$ is a shorthand for the predicted magnetic field value of
the $z$-component Gaussian process $\mathcal{GP}_{z}$ at $\bm{r}^{*}$.

The reason for the verbosity in Equation 3 is to allow an observation set for
optimizing hyperparameters
$D_{z}^{n_{2}}=(\bm{X}^{n_{2}},\tilde{\bm{Y}}_{z}^{n_{2}})$ and a separate
“inference set” $D_{z}^{n_{1}}=(\bm{X}^{n_{1}},\tilde{\bm{Y}}_{z}^{n_{1}})$
for computing the similarity weights when predicting the magnetic field at
some location $\bm{r}^{*}$. This is important because hyperparameter
optimization, typically done offline, scales with
$\mathcal{O}\left({n_{2}}^{3}\right)$ while inference scales with
$\mathcal{O}\left({n_{1}}^{2}\right)$. With this formulation, our map can
encode information from a large set of flights offline, then use a subset of
these observations during the actual inference to reduce computational load.
We leverage this idea in Section 4.4 to create and query our magnetic field
map.

###  3.2 Magnetometer Calibration

To calibrate our vehicle’s magnetometer, we use the model and iterative least
squares solver from Ref. Springmann and Cutler (2012) along with the two-step
calibration procedure from Ref. Wu et al. (2020). The magnetometer model from
Springmann and Cutler (2012) is repeated here

$\displaystyle\tilde{B}_{x}=\theta_{a}B_{x}+\theta_{x_{0}}+\eta_{x}$ (4)
$\displaystyle\tilde{B}_{y}=\theta_{b}(\tilde{B}_{y}\cos(\theta_{\rho})+B_{x}\sin(\theta_{\rho}))+\theta_{y_{0}}+\eta_{y}$
(5)
$\displaystyle\begin{split}\tilde{B}_{z}&=\theta_{c}(B_{x}\sin(\theta_{\lambda})+B_{y}\sin(\theta_{\phi})\cos(\theta_{\lambda})\\\
&+B_{z}\cos(\theta_{\phi})\cos(\theta_{\lambda}))+\theta_{z_{0}}+\eta_{z}\end{split}$
(6)

where ($\tilde{B}_{x}$, $\tilde{B}_{y}$, $\tilde{B}_{z}$) are measured
magnetic field values, ($B_{x}$, $B_{y}$, $B_{z}$) are the true magnetic
field, and the parameters to solve for are bias
($\theta_{x_{0}},\theta_{y_{0}},\theta_{z_{0}}$), scale factor
($\theta_{a},\theta_{b},\theta_{c}$), and non-orthogonality terms
($\theta_{\rho},\theta_{\lambda},\theta_{\phi}$).

Ultimately, we aim to find parameters
$\theta=[\theta_{a},\theta_{b},\theta_{c},\theta_{x_{0}},\theta_{y_{0}},\theta_{z_{0}},\theta_{\rho},\theta_{\lambda},\theta_{\phi}]^{\top}$
by minimizing the following cost function via an iterative non-linear least
squares solver

$\begin{split}\Delta B&=B_{R}^{2}-B^{2}\\\
&=B_{R}^{2}-(B_{x}^{2}+B_{y}^{2}+B_{z}^{2})\\\
&=B_{R}^{2}-g(\tilde{B}_{x},\tilde{B}_{y},\tilde{B}_{z},\alpha)\end{split}$
(7)

where $B_{R}$ is the reference magnetic field strength (taken from the World
Magnetic Model Team and Survey (2020)) and $g()$ is obtained by solving
Equations 4 \- 6 for ($B_{x},B_{y},B_{z})$.

The key difference from Springmann’s method Springmann and Cutler (2012) is to
estimate the nine $\theta$ parameters in two separate minimization steps like
Wu et. al. demonstrate in Wu et al. (2020). We found that this two-step
approach gave more consistent results for another magnetometer (not used in
this work) on our UAV. As such, we adopted the two-step calibration technique
for our primary magnetometer as well.

First, a simplified magnetometer model is created by setting all scaling terms
to 1 and non-orthogonalities to 0 leaving just the bias terms
($\theta_{x_{0}},\theta_{y_{0}},\theta_{z_{0}}$) in Equations 4 \- 6. The
first minimization is done on this simplified model to find an optimal set of
bias terms ($\theta_{x_{0}}^{*},\theta_{y_{0}}^{*},\theta_{z_{0}}^{*}$) which
are in turn used as initial conditions for the second optimization where
$\theta_{0}=$ $[1,1,1,$
$\theta_{x_{0}}^{*},\theta_{y_{0}}^{*},\theta_{z_{0}}^{*},$ $0,0,0]^{\top}$.
The nine optimal parameters from this second optimization are used as the
calibration terms for the magnetometer.

###  3.3 Performance Metrics

We use the root mean squared error (RMSE) of each GPR’s prediction
$\hat{m}_{x}$ against the corresponding component of the measured field
$\tilde{y}_{x}$ across a validation set. With this, the RMSE for
$\mathcal{GP}_{x}$ over some validation set
$D^{n_{v}}_{x}=(\bm{X}^{n_{v}},\tilde{\bm{Y}}_{x}^{n_{v}})$ is defined as

$RMSE_{x}=\sqrt{\frac{1}{n_{v}}\sum_{i=1}^{n_{v}}(\hat{m}^{\\{i\\}}_{x}-\tilde{y}^{\\{i\\}}_{x})^{2}}$
(8)

where $\hat{\bm{M}}_{x}^{n_{v}}=GPR_{x}(\bm{X}^{n_{v}})$ is the $x$ component
of the predicted magnetic vector, $\hat{m}^{\\{i\\}}_{x}$ is the $i$th
prediction in $\hat{\bm{M}}_{x}^{n_{v}}$, and $\tilde{y}^{\\{i\\}}_{x}$ is the
corresponding $i$th measurement in $\tilde{\bm{Y}}_{x}^{n_{v}}$.

Although we have three GPRs, it is sometimes convenient to summarize the
accuracy of their composite estimate as

$RMSE_{norm}=\sqrt{{RMSE_{x}}^{2}+{RMSE_{y}}^{2}+{RMSE_{z}}^{2}}.$ (9)

##  4 Experimental Procedure

This section introduces our quadrotor unmanned aerial vehicle (UAV), testing
locations, and flight trajectories used to construct and validate our magnetic
field maps.

###  4.1 Equipment, Facilities, and Setup

All tests were conducted at the Robot Fly Lab in the University of Michigan
Ford Motor Company Robotics Building. This indoor flight arena is equipped
with eight OptiTrack motion capture cameras that give a working volume of 4m
$\times$ 3m $\times$ 2.25m as shown in Figure 2. A ground station computer
connected to the OptiTrack system provides ground-truth position and attitude
estimates of the vehicle at 120Hz. The communication setup for the ground
station, pilot transmitter, and flight vehicle is the same as that depicted in
Figure 3 of Ref. Kuevor et al. (2021), but with one less BeagleBone Blue on
the UAV.

Figure 2: Flight vehicle Q1 in the Robot Fly Lab within the Ford Motor Company
Robotics Building at the University of Michigan.

The motion capture pose estimate is streamed to the flight vehicle in real-
time, but with a communication latency of about 40ms. For this, the vehicle’s
onboard estimate of roll and pitch is used for control while the remaining
states (position and heading) are taken from the 40ms-delayed motion capture
packets.

Time synchronization is important to properly associate each onboard
magnetometer measurement with a ground-truth pose when creating and validating
our magnetic field maps. During flights, our UAV’s BeagleBone Blue is
synchronized with the ground station laptop and ground station BeagleBone
Green using a Linux tool called chrony which was created by the Red Hat
Software company (https://chrony.tuxfamily.org/ accessed on February 2022). In
post-processing, this synchronization allows us to associate the non-delayed
motion capture data with the drone’s magnetometer readings. Additionally, we
use motion capture pose (position and attitude) estimates to rotate the
magnetometer data into the world frame and determine the location each
magnetic field measurement was collected.

This work uses time-invariant magnetometer calibration combining techniques
from Refs. Springmann and Cutler (2012) and Wu et al. (2020) (see Section
3.2). Magnetometer calibration is performed outdoors, just South of the
University of Michigan’s outdoor netted flight facility, M-Air. This location
is far enough away from any buildings that the magnetic field strength is
constant over a few meters and should be accurately reflected by the World
Magnetic Model (WMM) Team and Survey (2020). For the data gathered on
September 1st, 2022 (test series t6), the ambient magnetic field reference
term used in our calibration is $B_{R}=53.1351\mu T$ as taken from a WMM
online calculator for M-Air’s location at 42.294431∘N, 83.710442∘W, and 270m
above sea level.

Our flight vehicle “Q1” (Quadrotor one) is shown in Figure 3 and has an RM3100
magnetometer that is sampled at 200Hz and an MPU9250 IMU (gyroscope,
accelerometer, and magnetometer) sampled at 200Hz. This work made no use of
the magnetometer on the MPU9250 relying solely on the RM3100 for all magnetic
field measurements. The BeagleBone Blue (BBB) microprocessor handles
communication with the ground station via Xbee radio, logs all sensor data
locally, and commands the four motors to achieve the desired flight
trajectory. A Turnigy 4S 2650mAh 20C LiPo battery provides the main power on
Q1 while a Turnigy 2S 300mAh LiPo battery keeps the BBB powered on when
swapping 4S batteries between flights.

(a) Main components of Q1. $S_{2}$ = 2cm.

(b) Definition of $S_{1}$ and $S_{2}$. $S_{2}$ = 8cm.

Figure 3: Q1: The flight vehicle used for experiments in this work.

Figure 3(b) shows Q1 with different distances from the magnetometer to the
rest of the vehicle’s electronics. $S_{1}$ is defined as the distance from the
top of the battery mounting plate to the top of the BBB. This distance is
$S_{1}$ = 4.79cm on Q1 and will not change in any experiments through this
paper. $S_{2}$ is an adjustable distance defined from the top of the BBB to
the bottom of the magnetometer’s mounting plate that allows us to investigate
how the magnetometer’s proximity to other electronics affects the consistency
of magnetic field measurements on a quadrotor.

Finally, we placed a stationary RM3100 magnetometer in the corner of the
flight lab that gathers data every 10 seconds (0.1Hz). This stationary sensor
is on the ground at approximately (-3.5m, +2.5m) in the ($X_{W}$, $Y_{W}$)
frame defined in Figure 2 and was used to confirm that the ambient magnetic
field remained constant during our experiments.

Figure 4: Single-altitude lawnmower trajectory used in all trajectories in
this work.

###  4.2 Flight Profiles

The flight tests conducted for this work are all designed around single-
altitude scanning patterns that gather observations at a planar slice of the
working volume. Figure 4 shows the desired trajectory (blue, solid line) has
4m $x$-axis strides separated by 0.25m $y$-axis strides. Since all flights
begin with the drone at the origin (0m, 0m) of our working volume, the
trajectory in Figure 4 includes diagonal strides to leave and return to the
origin at the start and end of each flight test. Finally, some trajectories
gather observations at multiple altitudes in a single flight and include
$z$-axis strides between each planar trajectory. For multi-altitude tests, the
drone traverses a single diagonal stride from the origin to the corner at the
first altitude, and another diagonal stride from the corner to the origin at
the final altitude.

All linear strides are created from quintic spline trajectories that are hand-
tuned to enforce a 1.9m/s speed limit which is a compromise between the desire
for shorter flight tests (that allow for more frequent data collection) and
the limitations of our flight controller.

The 0.25m $y$-axis spacing is selected to enable a study on how the spatial
density of observations used to train a GPR-based map affects the accuracy of
the said map (Section 5.4). A similar study was done in ref. Akai and Ozaki
(2017) where magnetic field observations in a 3m $\times$ 3m $\times$ 2.2m
volume are hand-gathered at 0.2m increments in all axes. For our work, having
0.2m $y$-axis and 0.2m $z$-axis spacing that would scan all altitudes in our
working volume created a flight trajectory duration that exceeded the maximum
flight time of our vehicles. Thus, a 0.25m minimum separation distance is used
in this work instead.

####  4.2.1 Flight Profile Nomenclature

Throughout this paper, we refer to flights with ID tags like “tY_XX” where Y
refers to the flight test series and XX is the two-digit ID of the flight test
in that series. This allows the reader to reference our raw experimental data
to improve the reproducibility of our results. Appendix A has more details on
the matter.

###  4.3 Pre-Processing Magnetic Observations

When creating our hyperparameter ($D^{n_{2}}$) and inference sets
($D^{n_{1}}$), we first pre-process the magnetic field data gathered by the
RM3100.

Figure 5 shows data from the UAV’s RM3100 for a segment of time when the drone
is not moving and the motors are not spinning. The black dots are raw
magnetometer measurements that are mostly constant, but frequently have
spurious measurements that vary by $\pm 3\mu$T for each axis respectively. We
have seen such spurious measurements from three different RM3100 sensors; even
one mounted on a platform with no motors or ESCs (electronic speed
controllers). Other works using RM3100s have not reported such spurious
readings Hanley et al. (2021); Regoli et al. (2018). Thus, we believe this
could be a result of our 200Hz sampling rate (600Hz for the whole device to
achieve 200Hz per axis) and something with the RM3100’s firmware that does not
like such a fast sampling rate. Alternatively, this could be caused by some
poor voltage regulation on the BeagleBone Blue. We have not tried sampling
from other microprocessors to see if the BeagleBone is the cause of this
issue.

Nonetheless, all RM3100 data we use in this work is first put through a moving
median filter with a window size of 5. The resultant observations from the
median filter are shown in red in Figure 5. Clearly, there are still some
outliers in red, but increasing the window size of the moving median filter to
remove these outliers delays the signal enough to cause problems for state
estimation. Since we want to construct and validate maps with the same pre-
processing steps used when performing state estimation, we fix a window size
of 5 and work with the remaining outliers.

Figure 5: Data from a stationary RM3100 placed at the origin of our flight
arena and sampled at 200Hz. Spurious measurements are reduced with a moving
median filter with a window size of 5.

In addition to a median filter, we downsample the 200Hz of data gathered
during the flight test. Since the drone does not move very far in 1/200th of a
second, many observations are spatially redundant. As such, we temporally
downsample our observations to 2Hz, 4Hz or 10Hz when creating or validating
our GPR maps. Figure 4 shows the spatial distribution of 2Hz and 10Hz
downsampling during a flight test.

Since motion capture data is used to rotate the magnetometer data from the
body frame to the world frame, we perform one final check to ensure a timely
ground-truth pose estimate. Any temporally downsampled observation with a
motion capture pose older than 50ms is considered “stale” and replaced with
the nearest (in time) observation that has a more timely motion capture pose.
If the “fresh” replacement is already in the downsampled observation set, the
redundant observation is removed leaving one instance of it in the final set.

###  4.4 Creating and Querying the “Compromise” Map

As mentioned in Section 3.1, creating a GPR-based map requires each Gaussian
process $\\{\mathcal{GP}_{x}$, $\mathcal{GP}_{y}$, $\mathcal{GP}_{z}\\}$ to
have a set of hyperparameters optimized over an observation set
$\\{\bm{\Theta}^{*}_{x}(D_{x}^{n_{2}}),$ $\bm{\Theta}^{*}_{y}(D_{y}^{n_{2}}),$
$\bm{\Theta}^{*}_{z}(D_{z}^{n_{2}})\\}$ and an inference set
$\\{D_{x}^{n_{1}},D_{y}^{n_{1}},D_{z}^{n_{1}}\\}$ to compare the target
locations. This section explains how we use flight data to create the
hyperparameter observation sets (with $n_{2}$ observations) and then use
inference from an intermediate magnetic field map to make “compromise”
inference sets (with $n_{1}$ observations) where $n_{1}<n_{2}$. The goal of
the compromise map is to leverage the flight-by-flight variation in magnetic
field measurements (Section 5.2) by training hyperparameters on $n_{2}$
observations from many flights without incurring the computational cost of
having a large inference set size $n_{1}$.

The “hyperparameter observation set” for the $y$ component Gaussian process
$\mathcal{GP}_{y}$ is
$D_{y}^{n_{2}}=(\bm{X}^{n_{2}},\tilde{\bm{Y}}_{y}^{n_{2}})$ and similar for
$D_{x}^{n_{2}}$ and $D_{z}^{n_{2}}$. All three hyperparameter observation sets
share the same observed locations $\bm{X}^{n_{2}}\in\mathbb{R}^{3\times
n_{2}}$ taken from downsampled observations from any number of training
flights.

These observations are used to compute
$\\{\bm{\Theta}^{*}_{x}(D_{x}^{n_{2}}),\;\bm{\Theta}^{*}_{y}(D_{y}^{n_{2}}),\;\bm{\Theta}^{*}_{z}(D_{z}^{n_{2}})\\}$
by minimizing the negative log marginal likelihood of the observations over
the hyperparameters of each respective Gaussian process. The optimal
hyperparameters, along with their corresponding observation sets of size
$n_{2}$, will serve as an “intermediate” magnetic field map.

From here, we query our intermediate map at $n_{1}$ user-selected locations to
generate $n_{1}$ magnetic field observations for the compromise map. The
result of this step is a set of $n_{1}$ locations around the workspace that
make up $\bm{X}^{n_{1}}\in\mathbb{R}^{3\times n_{1}}$. The magnetic field
“measurements” that complete the inference set
$D_{x}^{n_{1}}=(\bm{X}^{n_{1}},\hat{\bm{M}}_{x}^{n_{1}})$ come from the
estimated magnetic field values at each location in $\bm{X}^{n_{1}}$ where

$\hat{\bm{M}}_{x}^{n_{1}}=\mathbb{E}\left[\mathcal{GP}_{x}\left(\bm{0},\;K_{\bm{\Theta}_{x}^{*}(D_{x}^{n_{2}})}(\bm{X}^{n_{2}},\bm{X}^{n_{1}})\right)\right]\in\mathbb{R}^{n_{1}}$
(10)

and similarly for the magnetic field “measurements” $\hat{\bm{M}}_{y}^{n_{1}}$
and $\hat{\bm{M}}_{z}^{n_{1}}$ for inference sets $D_{y}^{n_{1}}$ and
$D_{z}^{n_{1}}$ respectively. Here, the GPR maps estimate the magnetic field
at locations $\bm{X}^{n_{1}}$ by computing their similarity to locations
$\bm{X}^{n_{2}}$. The notation used in Equation 10 was introduced in Section
3.1.

At last, we have what we call our “compromise inference sets”
$(D_{x}^{n_{1}},D_{y}^{n_{1}},D_{z}^{n_{1}})$ (or just “inference sets”)
comprised of magnetic field estimates from the intermediate map at $n_{1}$
user-selected locations around our workspace.

Our construction of the compromise map gives two benefits. First, it allows
our GPR-based maps to avoid overfitting to a single training flight by
allowing large hyperparameter observation sets ($n_{2}$) that can incorporate
observations from several flights of a quadrotor (or potentially even flights
from multiple quadrotors for a multi-agent system). Next, it allows us to
perform inference in $\mathcal{O}({n_{1}}^{2}n_{q})$ rather than
$\mathcal{O}({n_{2}}^{2}n_{q})$ where $n_{q}$ is the number of points queried.

##  5 Results and Discussion

This section starts by explaining how the derivative gains from the PID
controller of our UAV’s autopilot create measurable magnetic noise (Section
5.1). Next, Section 5.2 introduces the unusual magnetic biases our UAV injects
into our measurements and how distancing the magnetometer from the electronics
improves the consistency of measurements. We then show that our compromise map
(which uses $n_{1}=511$ observations for inference) yields estimates within
0.013$\mu$T of our intermediate map (which uses $n_{2}=2001$ observations for
inference) in Section 5.3.

We follow with an analysis of the spatial density of observation points used
to train the GPR map in Section 5.4 and find that using observations within
0.55m of one another is sufficient to accurately represent the magnetic field
in our flight arena. Section 5.5 shows that it is equivalent to either create
a specialized map on the norm of the field or simply take the norm of
estimates from a vector-valued magnetic field map. Finally, Section 5.6
presents our consistency metric which is a tool to identify when a user can
rely on the predictions from their GPR-based magnetic field map.

###  5.1 Flight Controller Creating Magnetic Field Noise

Our vehicles use the rc_pilot_a2sys autopilot originally forked from the open-
source rc_pilot repository by James Strawson and Librobotcontrol
(https://github.com/StrawsonDesign/rc_pilot). rc_pilot_a2sys uses a four-stage
cascaded PID controller as explained in Section IV.B.1 of Romano et al.
(2022). Our work started with the gains used in Ref. Romano et al. (2022) but
adjusted them to reduce noise from the motors and ESCs.

Figure 6 shows magnetometer data from a segment of an experiment where the
flight vehicle (not Q1, but another vehicle of the same construction) was
commanded to a single-altitude scanning pattern (Figure 4) at 1.5m altitude.
The data in Figure 6 was rotated from the vehicle’s body frame to the world
frame using attitude estimates from motion capture cameras.

Figure 6(a) (t5_00) is the magnetometer data with the gains inherited from
ref. Romano et al. (2022). Here, the drone moves along the $x$ axis (4m
stride) of the flight space from 24s - 30s, and again from about 36s - 42s.
During these time segments, there is a large variance in the measured magnetic
field due to commands from the flight controller to the motors and ESCs.

Figure 6(b) (t5_01) shows another instance of the same trajectory, but with
new controller gains. The plots are temporally aligned to easily compare
corresponding flight segments. At a glance, it is clear that the variance in
the magnetic field measurements is much lower in Figure 6(b) than in Figure
6(a), but not completely gone. This is most evident in the “Mag Y - World”
plot of Figure 6(b) where the drone moves along the $x$ axis (4m stride) of
the flight space from 21s-26s and again from 34s - 40s.

(a) Test t5_00: Flight with “noisy” gains.

(b) Test t5_01: Flight with “quiet” gains.

Figure 6: Vehicle “Maggie” flying two instances of the same flight trajectory
showing how the flight controller can inadvertently create noise in the
measured magnetic field.

As explained in Section IV.B.1 of ref. Romano et al. (2022), the controller
used on this vehicle is a four-stage cascaded PID controller with an outer
loop that operates on position error and an inner-most loop controlling
angular rate error. The math for the inner loop is shown below in Equation 11

$\begin{split}\bm{\tau}_{\phi}&=K_{4,\phi}^{P}\,(\dot{\phi}_{d}-\;\tilde{\omega}_{x})\\\
&+K_{4,\phi}^{I}\,\int_{0}^{t}(\dot{\phi}_{d}-\;\tilde{\omega}_{x})\,ds\\\
&+K_{4,\phi}^{D}\,\frac{d}{dt}(\dot{\phi}_{d}-\;\tilde{\omega}_{x})\\\
\end{split},$ (11)

where $\bm{\tau}_{\phi}$ is the desired torque for roll $(\phi)$,
$\\{K_{4,\phi}^{P},K_{4,\phi}^{I},K_{4,\phi}^{D}\\}$ are the P, I, and D gains
for the fourth-stage controller along the roll axis, $\dot{\phi}_{d}$ is the
desired roll rate, and $\tilde{\omega}_{x}$ is the measured roll rate from the
gyroscope. There are similar fourth-stage PID loops for pitch ($\dot{\theta}$)
and yaw ($\dot{\psi}$) rates respectively.

The ‘noisy’ gains, inherited from Romano et al. (2022), used
$K_{4,\phi}^{D}=K_{4,\theta}^{D}=K_{4,\psi}^{D}=0.01$ while the modified
‘quiet’ gains set all three of these values to 0. This change alone is
responsible for the reduction in motor-induced noise in Figure 6.

The derivative term of the attitude rate controllers (e.g., Equation 11 for
roll rate) uses numerical derivatives of gyroscope measurements (Figure 7) as
part of their computation. As the quadrotor flies, the propellers induce
vibrations that make the gyroscope (and accelerometer) measurements rather
noisy. Taking numerical derivatives of these noisy gyroscope measurements
causes the motors to be commanded with high-frequency inputs that induce
measurable electromagnetic noise. We believe this problem is amplified as the
total amount of current pulled from the 4S batteries increases. This would
explain why the magnetometer variance is larger during the 4m, $x$-axis stride
(with a max commanded velocity of 0.75m/s for this set of tests) than during
the 0.25m, $y$-axis stride (max commanded velocity of 0.47m/s).

Figure 7: Raw gyroscope data from t5_00 . Vibrations from UAV propellers
create noisy gyroscope data. Differentiating this noisy data creates the high-
frequency magnetic noise seen in Figure 6(a).

Aside from changing controller gains, other possible solutions include
commanding less-aggressive maneuvers during flight tests, moving the
magnetometer further from the motors and ESCs (Section 5.2), or time-varying
magnetometer calibration Springmann and Cutler (2012). In addition, from some
preliminary data gathered when trying to understand this problem, we believe
that differences in motors and ESCs (even those of the same make and model)
may create different magnitudes of measurable magnetic noise.

The remainder of this paper uses only the ‘quiet’ controller gains. Note that
the analysis for this section used $S_{2}$ = 2cm (Figure 3(b)). In the
following section, we analyze how the distance of the magnetometer from the
motors and ESCs affects the variation in the measured magnetic field.

###  5.2 Varying $S_{2}$: Distance from Magnetometer to Electronics

Though the flight controller created motor-induced, high-frequency magnetic
field noise, there is another measurable magnetic anomaly at play with our
quadrotors that causes flight-by-flight variation in the measured magnetic
field. Here, we distinguish variance (the spread of a distribution of points)
from variation (more of an offset). We will reserve variance to describe the
spread of points around a probabilistic mean. For this, removing the
derivative gains from a particular PID loop of our flight controller reduced
the variance of motor-induced magnetic noise (Figure 6). To avoid ambiguity,
variation will be used to describe the differences in the magnetic field
between subsequent flights of the same trajectory. As we will show in this
section, the differences in measured magnetic field values are not obviously
spread around some average signal.

Recall that $S_{2}$ (Figure 3(b)) is the distance from the top of the
BeagleBone Blue (BBB) to the bottom of the magnetometer’s mounting plate. In
this section, we will vary the parameter $S_{2}$ to show that increasing the
distance of the magnetometer from the other electronics reduces the flight-by-
flight variations we see in the measured field. It is easier to explain these
anomalies graphically, which we will do in a moment, but first, we explain our
testing methodology.

####  5.2.1 Testing Methodology

Q1 was commanded to fly the same single-altitude scanning trajectory (Figure
4) at an altitude of 1.5m. Initially, in our investigation of these flight-by-
flight variations, we believed that each 4S LiPo battery could have a
different amount of “resting current” it provided to the motors to achieve the
same flight maneuver as another 4S battery. Though we could not conclusively
prove or disprove this hypothesis with our data, it did drive the design of
our experiments. For this, each battery started at full charge and Q1 flew the
same 1.5m scanning trajectory several times in a row. This allowed us to see
if the flight-by-flight variations were due to the battery’s voltage during a
flight. After a few repetitions with the first 4S battery, a new (fully
charged) 4S battery was used to power Q1 as it flew more repetitions of the
same trajectory. By testing different batteries, we could see if the cause of
the flight-by-flight anomalies was tied to which battery was powering the
motors and ESCs.

Three Turnigy 4S 2650mAh 20C batteries (denoted as #02, #04, and #14) were
used for this analysis. The charge/discharge and storage history of these
batteries are largely unknown as they have been used by multiple members in
the lab since they were purchased in May of 2021. The only thing we could
rigorously control for is fully charging each battery to 16.8V just before a
flight test.

For this analysis, each flight was downsampled to 10Hz or 2Hz to reduce
spatially redundant observations resulting in $\sim$900 or $\sim$180
observations respectively for each flight. Figure 4 shows the observation set
from one flight with 10Hz and 2Hz downsampling. Although 10Hz clearly shows
better coverage of the flight space, the RMSE values for this analysis changed
by at most 0.03$\mu$T on the $S_{2}$=8cm dataset when training on 10Hz vs 2Hz
downsampling sets. Thus, we sometimes use 10Hz downsampling to illustrate
specific points but typically use 2Hz downsampling for faster training and
validation of the magnetic field maps.

With this gathered data, we can train a GPR-based map of the magnetic field at
a 1.5m single-altitude slice of the working volume and compare the predictions
of that single-altitude map to the measurements gathered from the other
repetitions of the same trajectory. Since we assume the ambient magnetic field
is not changing throughout our experiments, we expect the magnetic field
measurements to be nearly identical between each flight.

####  5.2.2 Flight-by-flight variations

By training the map on a single repetition and validating on another, we get
data that looks like Figure 8 which has two types of plots: one with red and
blue lines with blue shading and another type of plot in grayscale.

For the first type (e.g., Figure 8(a)), each actual magnetometer measurement
(in the world frame) from the validation dataset is shown as a red cross. The
blue line is the GPR’s predicted mean at the location the drone was during
that timestamp in the validation flight. The blue shading around the solid
blue line depicts the two standard deviations ($2\sigma$) of the GPR’s
uncertainty with its prediction.

The second type of plot (e.g., Figure 8(b)) depicts a black dot as the error
between the red cross and blue line while the GPR’s $2\sigma$ uncertainty is
the gray shading. Here, we plot the absolute value of the error since the sign
of the error does not indicate anything of interest. The percentage in the
grayscale plots depicts how often the black dots (GPR error) are within the
gray shading (GPR’s $2\sigma$ uncertainty).

(a) $S_{2}$=2cm. Trained on t1_09 (#04). Validated on t1_05 (#14). Same data
as Figure 8(b)

(b) $S_{2}$=2cm. Trained on t1_09(#04). Validated on t1_05(#14). Same data as
Figure 8(a)

(c) $S_{2}$=2cm. Trained on t1_09(#04). Validated on t1_04(#14).

(d) $S_{2}$=2cm. Trained on t1_09(#04). Validated on t1_15(#02).

Figure 8: Flight-by-flight variations shown through three validation tests.

In Figure 8, the GPR map is trained on 2Hz-downsampled observations from a
single flight (t1_09) and validated against 10Hz-downsampled observations from
three flight tests from the same t1 flight series. These three validation
examples are generally representative of the different flight-by-flight
magnetic field anomalies we have seen. Figure 8(b) shows a relatively large
steady-state prediction error throughout, Figure 8(c) has a steady-state error
that changes suddenly partway through the flight, while Figure 8(d) is an
ideal case where the mapping and validation flights have agreeable
measurements.

This is why we refer to this anomaly as “flight-by-flight variation”. The type
of GPR error we see here across consecutive flights in the same test segment
(e.g., t1_XX) does not seem to be cleanly distributed around some mean bias
value. Instead, repetition t1_09 had similar magnetic observations as t1_15
(Figure 8(d)) yet starkly different observations than t1_05 (Figure 8(b)).
Further complicating the matter are cases like t1_04 (Figure 8(c)) with a
time-varying bias.

There is evidence in Figure 8 that suggests these magnetic anomalies are
caused by the ESCs and motors. The initial black dots in Figures 8(b) and 8(c)
have low errors before the anomalous bias takes effect. These beginning points
are observations gathered when the drone is at rest on the ground (before it
has taken off) and the bias sets in, if at all when the quadrotor is in the
air. This suggests that before the motors and propellers are spinning, the
observations from t1_05 and t1_04 are similar to those from t1_09.

This insight begs the obvious experiment of constraining the motion of the
quadrotor and sampling the magnetometer with and without the propellers
spinning. Unfortunately, we could not safely conduct such a test with our
current setup, but we plan to test this in the future.

Finally, it is possible that the ambient magnetic field in our flight arena is
changing over time. To check this, placed a stationary RM3100 in our workspace
for 3+ hours sampling data once every 10 seconds (0.1Hz). The raw data from
this experiment is shown in Figure 9. Here, we still have regular, spurious
measurements (like in Figure 5) but the main signal does not change by more
than 0.15$\mu T$ within each respective component of the measured magnetic
field. This is much smaller than the 1.8$\mu T$ of variation we see in the Z
component of Figure 8(b). Thus, we believe the flight-by-flight variations are
due primarily to UAV-induced noise and not by time-varying changes in the
ambient magnetic field.

Figure 9: Stationary RM3100 data sampled at 0.1Hz over 3+ hours during t6_XX
flight tests. $x$, $y$ or $z$ components of the magnetic field never change by
more than 0.15$\mu T$ respectively.

####  5.2.3 LiPo Batteries and Magnetic Variation

From the three cases in Figure 8, it is tempting to conclude that differences
in batteries are the cause of these flight-by-flight anomalies. Here, we see
that batteries #04 (used for t1_09) and #02 (used for t1_15) yield agreeable
observations. However, from Figures 8(b) and 8(c), it seems that battery #14
(used for t1_05 and t1_04) is disagreeable with the other two.

To investigate this idea, we flew 60 repetitions of the 1.5m-scanning
trajectory from Figure 4 over four test sessions (t1_XX, t2_XX, t3_XX, and
t4_XX) each corresponding to a different value of $S_{2}$ (2cm, 4cm, 6cm, and
8cm). For these tests, we used the same three batteries (#02, #04, and #14)
and flew consecutive repetitions of the scanning trajectory for each battery.
For each test session, we trained the GPR map on the first battery #04 flight
and validated the map on all flight tests of the same test session. Note that
the data in Figure 8 (tests t1_04, t1_05, t1_09, and t1_15) are four of the
$S_{2}$=2cm repetitions used in this larger analysis.

The results of this 60-flight analysis are summarized as a scatter plot in
Figure 10. The horizontal position of each point is meant to signify the value
of $S_{2}$ as 2cm, 4cm, 6cm, or 8cm. Any horizontal deviations from these
discrete values are only for visual clarity and do not reflect any deviation
in the actual value of $S_{2}$ when the data was gathered. The vertical
position is the vector RMSE of all three GPRs (Equation 9). Finally, the
colors and symbols identify which battery was used for each repetition with a
green $+$ for battery #04, an orange $\triangle$ for #02, and purple $\square$
for #14.

At $S_{2}$ = {2, 4, 6, 8}cm, each battery flew $N$ = {5, 4, 3, 8} repetitions
of the 1.5m scanning trajectory. For $S_{2}$ = 8cm, each battery started at
full charge, flew four consecutive reps, was re-charged, then flew four more
repetitions. We present the second charge of each battery in the $S_{2}$=8cm
dataset with a rotated $+$. $\triangle$, and $\square$ respectively.

Figure 10: $S_{2}$ versus vector GPR RMSE. At each $S_{2}$, the map is trained
on one flight from battery #04 (green +). For $S_{2}$={2,4,6,8}cm there are
$N$={5,4,3,8} reps per battery.

Figure 10 shows two things. First, increasing $S_{2}$ makes the observations
of each flight more consistent with one another. Since the GPR-based map (for
each $S_{2}$) is trained on a battery #04 flight, the scatter plot really
shows us the magnitude of the flight-by-flight variations relative to an
arbitrary battery #04 flight. As we increase $S_{2}$, the variations tend
towards a difference of about 0.3$\mu$T to 0.5$\mu$T.

Note that for each value $S_{2}$, there is a single green $+$ that is
validated on the same observations used to train the GPR map. Even so, we do
not see any of the validation tests going much below 0.3$\mu$T meaning this
may be something of a lower limit on how accurate our GPR maps can be with our
testing platform which would include inaccuracies due to sensor noise (a
function of our 200Hz sampling rate) and any electromagnetic noise the
quadrotor generates.

The second realization is that, unlike what the three examples in Figure 8 may
suggest, there is not a clear relationship between battery choice and flight-
by-flight variation. Although battery #14 (purple $\square$) tends to be most
disparate from battery #04 at $S_{2}$ of 2cm, 4cm, and 6cm, we see its
observations fall in line with the rest by $S_{2}$ = 8cm.

We take this moment to point out a flawed conclusion we stated in our previous
work. In Figure 10b of Kuevor et al. (2021), observations from a rectangular
trajectory (the first $\sim$500 points) are validated on a GPR-based map. In
that plot, there is about 2$\mu$T of an error on each of the $x$ and $z$ GPR
validation sets except before takeoff and after landing. In Kuevor et al.
(2021), we hypothesized that this error in $x$ and $z$ was caused by a
pitching misalignment in the magnetometer for that single flight (relative to
the pitch angle for all the other flights). However, if this were the case,
there likely would not be the brief moments of agreement before takeoff and
after landing as shown in Figure 10b of Kuevor et al. (2021). With what we
have learned in this work, we now believe this was actually an example of
these flight-by-flight variations.

For the work in ref. Kuevor et al. (2021), these anomalies may not have
happened as often due to us using a different instance of the M330 quadrotor.
Alternatively, it could be that the nature of attitude estimation (the focus
of Kuevor et al. (2021)) is fairly robust to 2-3$\mu$T variations between
flights if such anomalies do not significantly change the angle of the ambient
magnetic field (which has a magnitude of 40-70$\mu$T in our workspace). Thus,
it is likely that these variations were simply less salient in ref. Kuevor et
al. (2021) given our previous application of attitude estimation.

We believe that time-varying magnetometer calibration could address the
flight-by-flight variation issues presented here. In ref. Springmann and
Cutler (2012), they use measurements of the electric current near high-powered
devices to estimate and offset vehicle-induced magnetic fields from their
measurements. However, our vehicle does not have any electric current sensors,
so such a time-varying calibration solution is out of reach for our platform.
Instead, we address the problem with our “compromise” map.

###  5.3 Accuracy of Compromise Map

In this section, we propose a solution to the flight-by-flight variations that
incorporate data from several training flights into a single magnetic field
map. Since each flight has a chance of giving a different bias in the measured
magnetic field, the GPR-based map will overfit if trained on only a single
flight test. Thus, we instead, use $n_{2}$ observations from many flights that
span the working volume to create an “intermediate map” that optimizes
hyperparameters
$\\{\bm{\Theta}^{*}_{x}(D_{x}^{n_{2}}),\;\bm{\Theta}^{*}_{y}(D_{y}^{n_{2}}),\;\bm{\Theta}^{*}_{z}(D_{z}^{n_{2}})\\}$
and predicts with the $n_{2}$ observations. The compromise map then uses the
intermediate map’s estimates at $n_{1}$ user-selected locations to predict the
magnetic field using only $n_{1}$ points. The goal of this section is to see
how the accuracy of the compromise map varies as a function of $n_{1}$;
specifically, taking after ref. Akai and Ozaki (2017), to see how the spatial
density of the $n_{1}$ training points effects map accuracy.

####  5.3.1 Multi-Altitude Trajectories

This section uses data from eight different flight tests listed in Table 1
where the first four flights are typically used for training the
intermediate/compromise map while the remaining four are used for validation.
The trajectories listed in Table 1 are all multi-altitude tests where each
single-altitude slice is the same trajectory from Figure 4. “Lower Four alts.”
flies at altitudes $Z=[$-0.5, -0.75, -1.0, -1.25]m, “Upper Four alts.” at
$Z=[$-1.5, -1.75, -2.0, -2.25]m, “Scan-$\gamma$” at $Z=[$-0.5, -1.375,
-2.25]m, while “Scan-$\epsilon$” and “Scan-$\epsilon^{\perp}$” fly at
$Z=[$-0.75, -1.5, -2.0]m. The idea here is to gather redundant observations
throughout the working volume to help our GPR map learn the flight-by-flight
variations.

Additionally, Table 1 lists the number of 2Hz downsampled observations from
each flight. The number of 4Hz and 10Hz observations are approximately
$2\times$ and $5\times$ the values listed in Table 1.

Table 1: Separate flights are used for training and validation. All magnetic field observations were downsampled to 2Hz for training and validation. Number of observations each training flight contributes to the map is listed.  Flight | Flight | Training [2 Hz] | Validation [2 Hz]
---|---|---|---
Number | Description | (# Observations) | (# Observations)
t6_00 | Lower Four Alts. | 571 | –
t6_01 | Upper Four Alts. | 580 | –
t6_03 | Scan-$\gamma$ | 442 | –
t6_21 | Scan-$\epsilon^{\perp}$ | 408 | –
t6_04 | Lower Four Alts. | – | 573
t6_05 | Upper Four Alts. | – | 573
t6_06 | Scan-$\gamma$ | – | 441
t6_20 | Scan-$\epsilon$ | – | 384
Table 2: Multi-altitude map error. Trained on $n_{2}=408$ observations from flight test t6_21 (Scan-$\epsilon^{\perp}$). Flight ID and | GPR RMSE ($\mu$T) | Error within 2$\sigma$ (%)
---|---|---
Description | Norm|X|Y|Z | X|Y|Z
t6_04 (Lower Alts.) | 0.595 | 0.256 | 0.284 | 0.456 | 74.69 | 87.09 | 52.88
t6_05 (Upper Alts.) | 0.998 | 0.633 | 0.541 | 0.551 | 6.81 | 19.72 | 32.11
t6_06 (Scan-$\gamma$) | 0.642 | 0.362 | 0.318 | 0.424 | 56.01 | 86.62 | 68.25
t6_20 (Scan-$\epsilon$) | 0.556 | 0.229 | 0.293 | 0.413 | 56.25 | 47.40 | 39.06
Table 3: GPR error for multi-altitude intermediate map trained on $n_{2}=2001$ observations from tests t6_00, t6_01, t6_03, and t6_21. Flight ID and | GPR RMSE ($\mu$T) | Error within 2$\sigma$ (%)
---|---|---
Description | Norm|X|Y|Z | X|Y|Z
t6_04 (Lower Alts.) | 0.466 | 0.144 | 0.145 | 0.419 | 98.43 | 98.08 | 85.17
t6_05 (Upper Alts.) | 0.713 | 0.403 | 0.447 | 0.382 | 45.55 | 23.73 | 90.05
t6_06 (Scan-$\gamma$) | 0.378 | 0.184 | 0.155 | 0.292 | 97.96 | 96.83 | 96.83
t6_20 (Scan-$\epsilon$) | 0.406 | 0.111 | 0.222 | 0.321 | 99.48 | 92.97 | 93.75

####  5.3.2 Intermediate Map Accuracy

To start, we demonstrate the value of adding observations from many flights to
a GPR-based map. We call GPR maps that are trained on observations from
multiple flights an intermediate map.

Table 2 gives the performance of a multi-altitude magnetic field map trained
on $n_{2}=408$ (2Hz downsampling) observations from flight test t6_21 and
validated on observations from four different flight tests. The table is split
into three major columns that give the details of the validation flight and
the vector RMSE ($\mu$T) of all three GPRs along with the RMSE of each $x$,
$y$, and $z$ GPR respectively. The third major column quantifies how often (as
a percentage) an error data point lies within the 2$\sigma$ uncertainty of the
respective GPR. The performance metrics from the second and third major
columns are the same as those listed in the grayscale plots of Figure 8 (with
different training and validation flights in this section).

The underlying issue is the map will overfit to the observations gathered from
t6_21. In Section 5.2, we showed that there are flight-by-flight variations in
the measured magnetic field. By training a GPR-based magnetic field map on a
single flight, we prevent the map from learning to account for these flight-
by-flight variations. Note that training the map on $n_{2}=2038$ observations
(10Hz downsampling) from t6_21 gives norm RMSE values of 0.629$\mu$T,
1.034$\mu$T, 0.751$\mu$T, and 0.552$\mu$T when validating on t6_04, t6_05,
t6_06, and t6_20 respectively. Thus, simply increasing the number of training
samples from a single flight test does not improve the map’s performance.

By comparison, Table 3 trains the GPR hyperparameters (and performs inference)
on $n_{2}=2001$ observations (2Hz downsampling) from four different flight
tests: t6_00, t6_01, t6_03, and t6_21. Note, from Table 1, that we are now
training and validating on one instance each of all four types of flights
conducted for this analysis: lower four altitudes, upper four altitudes,
scan-$\gamma$, and scan-$\epsilon$. By comparing Table 2 to 3, we see that
adding observations from a variety of flights uniformly reduces RMSE and
increases the frequency that error falls within each GPR’s 2$\sigma$ error
(usually caused by an increase in each GPR’s uncertainty).

Additionally, training on $n_{2}=9999$ observations (10Hz downsampling) from
the four training flights gives norm RMSE values of 0.497$\mu$T, 0.686$\mu$T,
0.406$\mu$T, and 0.408$\mu$T when validating on t6_04, t6_05, t6_06, and t6_20
respectively. Again, a bit to our surprise in this second case, adding more
observations from the same ensemble of training flights does not necessarily
improve the map’s performance. However, it is clear that training on multiple
flights (Table 3) is better than training on a single flight (Table 2).

Aside from the 2Hz vs. 10Hz comparisons done in this section, we do not aim to
directly address ways to reduce the cost of training hyperparameters. Instead,
we focus on the inference cost $\mathcal{O}\left({n_{2}}^{2}\right)$ with what
we call a compromise map.

The idea here is to query the intermediate map (which uses $n_{2}$
observations for inference) at $n_{1}$ user-selected locations
$\bm{X}^{n_{1}}\in\mathbb{R}^{3\times n_{1}}$ in the working volume. The
intermediate map’s estimates of the magnetic field at these $n_{1}$ locations
give us a set of “measurements” $\hat{\bm{M}}^{n_{1}}\in\mathbb{R}^{3\times
n_{1}}$ which (together with $\bm{X}^{n_{1}}$) form the “observations” used by
the compromise map to perform inference. Note that the compromise map will use
the same hyperparameters from the intermediate map for its inference. Section
4.4 gives a more detailed explanation of the process of creating a compromise
map from an intermediate map.

####  5.3.3 Intermediate Map vs. Compromise Map

We now compare the accuracy of the intermediate map, which uses $n_{2}=2001$
observations (2Hz downsampling) to perform inference, to that of the
“compromise” map which uses only $n_{1}=511$ inference points. Recall that
both use the same sets of hyperparameters optimized over $n_{2}=2001$
observations.

The $n_{1}=511$ user-selected locations for this analysis are chosen as
follows. Points are distributed evenly through the [-2, 2]m $x$ axis span,
[-1.5, 1.5]m $y$ axis span, and [-2.25, -0.5]m $z$ axis span of the working
volume. A compromise map location is selected every 0.5m, 0.5m, and 0.25m for
the $x$, $y$, and $z$ axes respectively. In total, this gives 504 locations
within the working volume. The remaining 7 points are evenly spaced from the
ground to an altitude of 0.5m so the compromise map has some observations
during the takeoff and landing sequence (which, for all our flights, are above
the origin). The number $n_{1}=511$ is an important constraint for another
toolbox we used for position localization. Thus, this spatial discretization
(0.5m $\times$ 0.5m $\times$ 0.25m) became a common test state in our work.

We can see these trends better in Figure 11 which is similar in style to the
grayscale plots from Figure 8, but with the subplots as three columns rather
than as three rows. This format fits more figures on a single page to more
easily compare the intermediate and compromise maps.

Figure 11 gives plots for the intermediate map on the left column and the
compromise map on the right. Here, we see the compromise map has both
quantitatively and qualitatively similar performance to the intermediate map
despite using nearly a fourth of the points for inference. By comparing the
norm RMSE values (in the titles) across the two columns of Figure 11, we see
the two maps are never off by more than 0.013$\mu$T (13 nanoTesla) in norm
RMSE. This difference is within the noise of the magnetic field measurements
of our platform.

For reference, training a GPR map on 408 observations (2Hz downsample) from
t6_21 and validating this map on the same 408 observations gives RMSE values
of (0.196$\mu$T, 0.089$\mu$T, 0.108$\mu$T, 0.136$\mu$T) for norm, $x$, $y$,
and $z$ RMSE values respectively. This means that on our flight platform, even
when the map has overfit for the exact observations it will be validated on,
it still cannot discern differences of 0.013$\mu$T.

Thus, we take the differences between the intermediate map ($n_{2}=2001$) and
the compromise map ($n_{1}=511$) in Figure 11 to be negligible and assert that
their RMSE performance is effectively identical.

Of course, a lot of valuable information can be lost by simply comparing RMSE
values. However, visual comparison across the two columns of Figure 11 further
emphasizes that both the intermediate and compromise maps yield very similar
prediction results. The key difference is in the uncertainty of the two maps.

It’s easiest to see this in row three (Figures 11(d) vs. 11(e)) where the
compromise map (right) has spikier gray shading (2$\sigma$ uncertainty) than
the intermediate map (left). The compromise map’s increased uncertainty is
primarily due to our choice of prediction locations
$\bm{X}^{n_{1}}\in\mathbb{R}^{3\times n_{1}}$ which have 504 points selected
from a (0.5m $\times$ 0.5m $\times$ 0.25m) spatial discretization of our
working volume. Recall that each altitude of our validation trajectories
traverses lanes separated by 0.25m (Figure 4). Thus, anytime the compromise
map is queried at a location between the (0.5m $\times$ 0.5m) $x$-$y$ points,
it will report a higher uncertainty in its predicted magnetic field estimate.

The overall increased uncertainty of the compromise map also explains why the
respective GPR errors more frequently fall within two standard deviations of
their respective uncertainties. Said another way: the plots on the right
column of Figure 11 will almost always have higher red percentages than the
comparative GPR on the left column given our selected spatial discretization
of (0.5m $\times$ 0.5m $\times$ 0.25m).

In general, we expect the overall uncertainty of the compromise map to
decrease with increased spatial density in our selection of locations for
$\bm{X}^{n_{1}}$. The next section will further investigate how the spatial
density of the points $\bm{X}^{n_{1}}$ affect the performance of the
compromise map.

(a) $n_{2}$ for t6_04. No outliers.

(b) $n_{1}$ for t6_04. No outliers.

(c) $n_{2}$ for t6_05. Outliers {[2.89], [], [2.42]}$\mu$T.

(d) $n_{1}$ for t6_05. Outliers {[2.49],[],[2.39]}$\mu$T.

(e) $n_{2}$ for t6_06. Outliers {[2.41],[],[]}$\mu$T.

(f) $n_{1}$ for t6_06. Outliers {[2.36], [], []}$\mu$T.

(g) $n_{2}$ for t6_20. No outliers.

(h) $n_{1}$ for t6_20. No outliers.

Figure 11: Performance of the intermediate map (left column) and compromise
map (right column). Outliers larger than 2$\mu$T are listed for {
[$\mathcal{GPR}_{x}$], [$\mathcal{GPR}_{y}$], [$\mathcal{GPR}_{z}$] } in each
subplot. Empty brackets indicate no hidden outliers for the respective
$\mathcal{GPR}$.

###  5.4 Compromise Map - Spatial Density Analysis

This section seeks to understand how the accuracy of the compromise map
changes as a function of the spatial density of the locations used in
$\bm{X}_{1}$. Generally, we expect GPR error to increase, the number of error
points captured by 2$\sigma$ uncertainty to increase, and the computation time
to decrease as the training set points become more sparse. In this section, we
will only analyze the first of these (GPR accuracy).

In the last section, our $n_{1}=511$ compromise training points came from a
custom spatial density of (0.5m $\times$ 0.5m $\times$ 0.25m) which gave 504
locations, plus an additional 7 locations to have observations along the
takeoff and landing segment for each flight. In this section, we will work
with a uniform spatial density in all directions ($S$ $\times$ $S$ $\times$
$S$). Additionally, to simplify matters a bit, we will use a naive linear
spacing of [$min:S:max$] for the locations along each respective spatial axis
of the lab. This means that different spatial densities $S$ can yield the same
number of points $n_{1}$ but at different locations. Recall that our working
volume has spatial limits of [-2, 2]m in $x$, [-1.5, 1.5]m in $y$, and [-2.25,
-0.5]m in $z$.

In this study, we vary $S$ from 0.2m to 1m to emulate the study done in ref.
Akai and Ozaki (2017). However, our study will quantify the RMSE of the
magnetic field map rather than relying exclusively on visual comparisons of
the map. For the 17 values of $S$ tested (in increasing order), $n_{1}$ =
[3031, 1775, 931, 655, 447, 259, 259, 199, 133, 112, 97, 97, 79, 67, 47, 47,
47]. Recall, seven of the $n_{1}$ values listed here are for the takeoff and
landing sequence.

Figure 12 shows the norm RMSE of the compromise map (validated on t6_04,
t6_05, t6_06, t6_20) as a function of the spatial density term $S$. A zoomed
version of the initial segment of data is embedded in the same figure. Here we
see that the norm RMSE is fairly constant for values $S\leq 0.5$m with a small
spike in norm RMSE at $S=0.45$m. This spike is caused by our naive linear
spacing which occasionally causes poor sampling along altitudes. At $S=0.45$m,
we get altitudes of $z=$ [-2.25, -1.8, -1.35, -0.9]m. Recall from Section
5.3.1 that some of our validation flights have measurements as low as $z=$
-0.5m.

Figure 12: GPR norm RMSE fairly constant for $S\leq 0.55$m. Flights t6_04 (red
$\square$) and t6_06 (black $\times$) have observations at $z=-0.5$m causing
higher errors for certain values $S$.

Combining these sources of information, we can see that the black-$\times$
(t6_06; Scan-$\gamma$) and red-$\square$ (t6_04; Lower four alts) suffer the
worst increases at $S=0.45$m. Since the compromise map’s observations only go
as low as $z=-0.9$m, both these trajectories have higher GPR RMSE when
validated at their $z=-0.5$m observations. Meanwhile, blue-$\circ$ (t6_20;
Scan-$\epsilon$) shows little change at $S=0.45$m since all its observations
are away from the extremes of the working volume. The other transient spikes
in Figure 12 (e.g., $S=0.6$m) are caused by similar $z$-axis sampling from our
naive linear spacing. t6_04 and t6_06 continue to be the most sensitive to
certain values $S$ when the map has no observations near $z=-0.5$m.

Despite these transient spikes, there are two clear trends we can see from
Figure 12. First, for values of $S\leq 0.55$m (with $S=0.45$m as an
exception), the norm RMSE is relatively insensitive to changes in $S$. This
agrees quite well with the ref. Akai and Ozaki (2017) where they show
(visually) that their magnetic field map is qualitatively similar for values
of $S=$ [0.2, 0.4, 0.6]m. In Akai and Ozaki (2017), their analysis led them to
use 0.6m as a standard distance between observations in another experiment
they conduct.

For comparison, we briefly summarize the spatial density of observations used
by some other indoor magnetic field mapping works. Some works map along a
single axis ($p=1$) such as Haverinen and Kemppainen (2009) that separates
observations by 0.04m and 0.1m on a ground robot and chest-mounted pedestrian
setup respectively through hallway networks of up to 350m in length.
Additionally, Suksakulchai et al. (2000) makes a $p=1$ map with 0.2m between
each observation and Akai and Ozaki (2015) make a $p=1$ map with 0.1m between
observations.

For $p=2$, planar maps Li et al. (2012) create a constant-altitude planar map
using $S=0.305$m separation in $x$ and $y$ for their coarse grid and
$S=0.005$m for fine grid spacing. This mapping is done in a 2.1m$\times$2.1m
room. Ref. Vallivaara et al. (2011) uses a type of occupancy grid for their
$p=2$ magnetic field SLAM solution with grid size 0.05m$\times$0.05m. The
solution ignores all previous magnetic field measurements further than 0.5m
from the vehicle’s current position. Ref. Wu et al. (2019) makes $p=2$ maps
with lanes separated by 0.38m. Ref. Hanley et al. (2021) creates $p=2$ maps at
different heights along hallways in three separate buildings. The authors use
a longitudinal separation distance of about 0.4m (varies per building) and a
vertical separation of 0.08m between their observations.

Finally, Ref. Brzozowski et al. (2016) makes a $p=3$ “map” of a hallway (no
explicit interpolation of observations) with (0.33m $\times$ 0.45m $\times$
$S^{?}$) where the vertical spacing $S^{?}$ is not explicitly listed in their
work, and Lipovský et al. (2021) makes $p=3$ maps with $S=0.5$m uniform
spacing in all axes. Our previous work used planar lawnmower patterns like in
Figure 4, but with 0.5m spacing in the $y$ axis (instead of the 0.25m spacing
used in this work) Kuevor et al. (2021). Furthermore, we used a vertical
spacing of 0.75m in our training sets to investigate the accuracy of $p=3$
magnetic field maps when interpolating between and extrapolating outside of
pre-mapped altitudes.

In short, for 1D, 2D, or 3D indoor magnetic field maps, we have not seen any
previous works using separation distances $S$ larger than 0.6m. The exceptions
are Akai and Ozaki (2017) where the authors explicitly study a reasonable
upper limit on $S$ and our work Kuevor et al. (2021) which studies the
interpolative and extrapolative limits of $p=3$ GPR maps.

We believe the empirical upper bound of $S=0.6$m described above might be
driven by proximity to walls and the floor. If this is the case, $S\leq 0.6$m
might only be the upper bound for indoor regions that are within a couple
meters from some wall or the floor. We acknowledge that this includes nearly
every room or hallway in most buildings, but a large, multi-story atrium might
be an exception. Portions of such a large, open, yet indoor volume may have
lower spatial variation and can be accurately mapped with values $S>0.6$m. The
size of the two works that explicitly study spatial density indoor magnetic
field mapping (3m $\times$3m$\times$2.2m for Akai and Ozaki (2017) and
4m$\times$3m$\times$2.25m in this work) might be too small (i.e., too close to
the walls and floor) to test this idea.

It is also possible that the floor, and not walls, are the dominating factor
for spatial variation in indoor magnetic fields. In Hanley et al. (2021),
Hanley et al. show that that the magnetic field near the floor (within 0.5m
from the ground) of three university buildings is most distinct from the field
at other altitudes up to 2m. This might explain why the other works in our
spatial density summary, which mostly use ground robots and pedestrians that
remain “near” the floor, have $S\leq 0.6$m.

The second conclusion, the complement of the first, is that values of $S\geq
0.6$m start to show steady increases in norm RMSE. It is not clear what trends
hold for $S>1$m, but such separation distances are somewhat degenerate given
the size of our working volume. Further, it is unlikely for the norm RMSE
values to return to their $S=0.2$m for values of $S>1$m (which is what we are
most interested in).

This analysis is what led us to use $S=0.5$m for the $x$-$y$ separation
distance for the $n_{1}=511$ compromise map used elsewhere in this paper. We
chose $S=0.25$m for $z$ because we found (empirically) that the GPR’s
performance is most sensitive to omissions in observed altitudes. Recent works
on indoor magnetic field mapping have shown similar trends: that indoor
magnetic fields change noticeably as a function of altitude Kuevor et al.
(2021); Hanley et al. (2021). For our analysis, similarly for Kuevor et al.
(2021); Hanley et al. (2021), it is not clear how much this altitude-based
sensitivity is driven by the fact that all our flights (training and
validation) are comprised of many single-altitude slices.

This section has shown that our compromise map can perform accurate magnetic
field inference with fewer prediction points $n_{1}<n_{2}$. Remember the goal
here was to allow the GPR-based magnetic field map to learn the flight-by-
flight variations across several training flights without incurring high
computation costs when performing inference. However, some validation flights
like t6_05 (Figures 11(c) and 11(d)) still have steady-state errors that even
the compromise map (or the intermediate map) is unable to resolve given a set
of training flights.

In the following section, we compare two methods of mapping the norm of the
magnetic field in our workspace to investigate if composing the three outputs
from a $3\rightarrow 3$ map gives a less-accurate norm estimate than training
a specialized $3\rightarrow 1$ map on the magnetic field norm.

###  5.5 Creating a $3\rightarrow 1$ map from estimates of a $3\rightarrow 3$
map

This section compares two ways of estimating the norm of the magnetic field in
a workspace. The first takes the three estimates from our $3\rightarrow 3$ map
and computes the norm of these estimates
$\sqrt{\hat{m}_{x}^{2}+\hat{m}_{y}^{2}+\hat{m}_{z}^{2}}$. The alternative is
to create a new $3\rightarrow 1$ map, called $GPR_{nrm}$, trained on the norm
of our measurements

$\begin{split}\tilde{y}_{nrm}&=\sqrt{\tilde{y}_{x}^{2}+\tilde{y}_{y}^{2}+\tilde{y}_{z}^{2}}\in\mathbb{R}\\\
\end{split}$ (12)

to create hyperparameter ($D_{nrm}^{n_{2}}$) and inference ($D_{nrm}^{n_{1}}$)
sets. The scalar predictions from this new GPR are defined as
$\hat{m}_{nrm}\in\mathbb{R}$.

The goal here is to see if a specialized $3\rightarrow 1$ map can more
accurately estimate the magnetic field norm than a composition of the
estimates from a $3\rightarrow 3$ map. This is an important distinction since
some works create maps on special components of the ambient magnetic field
like a $2\rightarrow 2$ map on the horizontal and vertical components of the
field Le Grand and Thrun (2012); Vandermeulen et al. (2013) or a $1\rightarrow
1$ map on the heading (or declination) of the field Suksakulchai et al.
(2000). Here, we aim to check if composing the three estimates from a
$3\rightarrow 3$ map accumulates error in estimating the norm of the field in
ways that a specialized $3\rightarrow 1$ map, which estimates the norm
directly, would not.

The procedure for training hyperparameters and creating the $n_{1}=511$
compromise map for our new $GPR_{nrm}$ is the same as described in Section
4.4, but with measurements as defined in Equation 12.

The RMSE accuracy of our maps requires a prediction from a GPR (or a set of
GPRs) at location $\bm{r}^{*}$ and a measurement taken at the same location.
We refer to the prediction given by our typical $3\rightarrow 3$ map as our
vector prediction

$\begin{split}\hat{\bm{m}}_{vec}(\bm{r}^{*})=\begin{bmatrix}\hat{m}_{x}(\bm{r}^{*})\\\
\hat{m}_{y}(\bm{r}^{*})\\\
\hat{m}_{z}(\bm{r}^{*})\end{bmatrix}&=\begin{bmatrix}\mathcal{GPR}_{x}(\bm{r}^{*})\\\
\mathcal{GPR}_{y}(\bm{r}^{*})\\\
\mathcal{GPR}_{z}(\bm{r}^{*})\end{bmatrix}.\\\ \end{split}$ (13)

Now we define two error terms meant to compare the accuracy of estimating the
magnetic field norm in our workspace. One is for the accuracy of our
$3\rightarrow 3$ vector map ($vec$) while the other is for the $3\rightarrow
1$ norm map ($nrm$).

$\begin{split}e_{vec\\_nrm}&=\sqrt{\hat{m}_{x}^{2}+\hat{m}_{y}^{2}+\hat{m}_{z}^{2}}-\tilde{y}_{nrm}\\\
e_{nrm\\_nrm}&=\hat{m}_{nrm}-\tilde{y}_{nrm}.\\\ \end{split}$ (14)

The error terms in Equation 14 are used to create RMSE metrics (Equation 8)
for the accuracy of each respective GPR against some validation set.

Table 4: Norm ($nrm$) GPR versus vector ($vec$) GPR on estimating norm of magnetic field. Flight ID and |  RMSE ($\mu$T)
---|---
Description | $e_{vec\\_nrm}$ | $e_{hrz\\_nrm}$
t6_04 (Lower Four Alts.) | 0.391 | 0.390
t6_05 (Upper Four Alts.) | 0.256 | 0.257
t6_06 (Scan-$\gamma$) | 0.245 | 0.246
t6_20 (Scan-$\epsilon$) | 0.301 | 0.301

Table 4 computes RMSE on the error metrics of Equation 14 for the four
validation flights we used in the previous section. From this, we see that
both methods of estimating the norm of the magnetic field are equivalent.
There is never more than 1nT of difference between the two RMSE values which
is well within the noise floor of the PNI RM3100 (even without accounting for
noise generated by the drone).

Our result is reassuring in that there is no need to train a fourth GPR
$GPR_{nrm}$ if a user wants the three vector components and the norm of the
magnetic field in their workspace. It is not immediately clear if this result
extends to angle-based compositions of the magnetic field map like declination
Suksakulchai et al. (2000) or inclination of the ambient field.

Our last tool in this paper, which we call the “consistency metric”, uses the
frequency of error falling with a GPR’s 2$\sigma$ interval to indicate if a
validation flight is consistent with a given GPR map.

###  5.6 Consistency Check

Despite our efforts to reduce magnetic noise from motors and ESCs (Sections
5.1 and 5.2) and incorporate flight-by-flight variations into a compromise map
(Section 5.3), there may still be flight tests with significantly different
magnetic field observations than what a given GPR map can predict. Such cases
can arise from the flight-by-flight variations of the same vehicle (e.g.,
Figures 11(c) and 11(d)), but may also come into play if different platforms
are used for mapping and utilization (e.g., GPR map trained on drone data but
utilized for indoor pedestrian localization). As such, it is important for
users to know when they can rely on the predictions from their GPR-based
magnetic field map, and when it may be better to leverage a different sensing
modality for their state estimates.

For this, we introduce the concept of “consistency” between a GPR-based map
and a validation flight. Given a $p\rightarrow m$ magnetic field map, we say
that a validation flight is consistent with the map if all $m$ GPRs
respectively capture 96% of their error within two standard deviations
(2$\sigma$) of their uncertainty.

Of course, the 96% threshold is inspired by the expected number of points that
fall within two standard deviations of a univariate normal (Gaussian)
distribution. Additionally, we note that a 96% threshold works well for our
test platform, our choice of GPR kernel (squared exponential), the way we
optimize hyperparameters (minimization of log marginal likelihood; Section
3.1), and our target application of position localization. Although these
kernel and optimization methods are common, we have not investigated how GPR-
based maps generated with different kernels or different hyperparameters
behave in regard to our notion of consistency.

To demonstrate, we create a $n_{1}=511$ compromise map of 2Hz-downsampled
observations from the following seven flight tests: t6_0, t6_1, t6_3, t6_4,
t6_5, t6_6, t6_16. We validate this compromise map on 10Hz samples from four
flight tests (t6_09, t6_11, t6_15, t6_18) to illustrate how we use the
consistency metric. Note that we are training and validating on different
t6_XX tests compared to the previous sections. The training flights used here
are from an analysis of magnetic field position localization which we will
explore in a future paper.

(a) t6_09. No hidden outliers.

(b) t6_11. No hidden outliers.

(c) t6_15. Outliers {[2.16, 2.56], [2.10], []}$\mu$T.

(d) t6_18. No hidden outliers.

Figure 13: Consistency metric (red percentage) across four validation flights.
Some error points are larger than 2$\mu$T and are cut off in the graphs. Such
outliers are listed for { [$\mathcal{GPR}_{x}$], [$\mathcal{GPR}_{y}$],
[$\mathcal{GPR}_{z}$] } in the respective subplot caption. Empty brackets
indicate no hidden outliers for the respective $\mathcal{GPR}$.

Figure 13 shows the GPR’s performance on the four selected validation flights.
We start with a case that easily passes our consistency test (Figure 13(a)).
Test t6_11 (Figure 13(b)) fails the consistency test since the
$\mathcal{GP}_{z}$ does not capture enough of its error within its $2\sigma$
uncertainty. Note that even though $x$ and $y$ GPRs have a great understanding
of their error, the poor performance of the $\mathcal{GP}_{z}$ alone means
that test t6_11 is inconsistent with this $3\rightarrow 3$ map. However, if
the user needed only the outputs of the $\mathcal{GP}_{x}$ and
$\mathcal{GP}_{y}$ (say for localization using magnetic field heading
Suksakulchai et al. (2000)), then the poor performance of $\mathcal{GP}_{z}$
may be irrelevant and the map could still be useful for such an application.

Figure 13(c) shows an example of the flight-by-flight variation changing
partway through a validation flight. In this case, the first 1220 validation
points easily pass our consistency metric for all three GPRs. From 1220 to
1335, things get worse and $\mathcal{GP}_{x}$ barely fails our consistency
test. After data point 1335, all three GPRs fail the consistency check.

This example allows us to emphasize an important point about how we use the
consistency metric. Since most of test t6_15 (up to point 1220) is consistent
with the compromise map, computing our 2$\sigma$ percentage across the entire
validation set gives (97%, 95%, 95%) for each GPR respectively. Working with
these numbers alone (without having Figure 13(c) for context), it’s tempting
to say all of t6_15 is consistent with the map (although barely). The reality
is that since the flight-by-flight variations can change bias values partway
through a single flight test, it is sometimes necessary to see if a portion of
a flight test fails the consistency test rather than checking the 2$\sigma$
metrics across the whole flight.

Finally, validating on t6_18 gives Figure 13(d) which is another gray area in
our use of the consistency check. Here, the first 1242 validation points
undoubtedly pass the consistency test. However, the remaining points fail the
consistency check on the z-component of the magnetic field. Note that the
$\mathcal{GP}_{y}$ plot in Figure 13(d) also changes behavior near index 1380,
but our analysis focuses on $\mathcal{GP}_{z}$ since it fails the consistency
check before $\mathcal{GP}_{y}$ does.

96% should be taken as a “rule of thumb” rather than a magic threshold where
something fundamentally different occurs with our construction of magnetic
field maps. With this in mind, it is important to consider how the map will be
used when deciding whether a validation test should be rejected on the basis
of the consistency check alone. Here, for t6_18 in Figure 13(d), we also
include the RMSE values (in blue) for all points before and after 1242 to
illustrate the accuracy of our compromise map before and after the mid-flight
bias shift.

For our analysis in Ref. Kuevor et al. (2021), though we did not realize it at
the time, these flight-by-flight variations were still present. Yet, we were
still able to use our magnetic field maps (at the time, trained on a single
flight test rather than our new compromise map method) to perform good
attitude estimates. We believe this is due to the nature of attitude
estimation which compares the angle of a measured vector against its
corresponding reference vector. In this case, the flight-by-flight variations
did not significantly change the angle of the magnetic field reference vectors
given by our GPR maps.

However, our preliminary results on magnetic field position localization show
that inconsistencies like those at the end of t6_15 (Figure 13(c)) and the end
of t6_18 (Figure 13(d)) is enough to confuse a particle filter and cause
noticeably larger estimation errors. The use of GPR-based magnetic field maps
for position localization (via a particle filter) has been the motivation for
much of this paper. However, if the intended application is more robust to
“small” inaccuracies of the GPR map, then a more lenient consistency check
threshold may work.

##  6 Conclusions and Future Work

This work presents practical methods to address vehicle-induced noise when
creating and validating GPR-based magnetic field maps with a UAV. First, we
showed that certain portions of a vehicle’s PID controller can increase the
measurable magnetic noise when commanding the motors and ESCs on a flight
vehicle.

Next, we introduced the flight-by-flight magnetic field variations present on
our quadrotor and presented two ways to address them. The first involves
moving the magnetometer further from the electronics on the vehicle while the
second aims to teach the GPR-based map about the flight-by-flight variations
through a “compromise” map. We show that our compromise map has similar
accuracy to the “intermediate” map which uses four times as many reference
points to estimate the magnetic field through in our working volume. We found
that a compromise map trained on observations spaced no more than $S=0.55$m
give comparable accuracy.

Finally, we note that our proposed methods of reduction and incorporation do
not encompass all possible mapping errors. As such, we introduce a
“consistency check” as a rule of thumb to help experimentalists quickly assess
if magnetic field observations from a particular test dataset are consistent
with their magnetic field map.

In the future, we aim to use our new mapping techniques to demonstrate the 3D
position localization of a UAV using indoor magnetic field maps. Per Section
5.5, we aim to conduct studies similar to Le Grand and Thrun (2012); Haverinen
and Kemppainen (2009); Vandermeulen et al. (2013) that compare the position
localization accuracy of magnetic maps with different output dimensions $m$.
The proposed work would complement our previous work on UAV attitude
estimation using indoor field maps Kuevor et al. (2021) which we will consider
revisiting with our new mapping methods. Additionally, investigating time-
varying magnetic fields (due to the position of large, metallic objects like
elevators and doors) is important and can be complicated by vehicle-induced
magnetic variations that may be confused for a change in the ambient magnetic
field. So far, our work assumed the indoor magnetic field is constant. This
tends to be a good assumption unless large metal objects (like doors and
elevators) are moved in the mapped area Haverinen and Kemppainen (2009); Solin
et al. (2015); Michel et al. (2018). As such, we aim to extend our methods
into time-varying magnetic field mapping that uses our consistency check
(Section 5.6) to decide when a portion of a map may be out of date and needs
to be updated.

Finally, a component-wise causal analysis of observed flight-by-flight
variations may better help us understand measurable bias and noise sources to
modify or eliminate entirely. For example, it is possible that simply using
different brands/models of ESCs and motors could make the UAV more
magnetically quiet.

Conceptualization, all; methodology, P.E.K., E.M.A. and J.W.C.; software,
P.E.K.; validation, P.E.K.; formal analysis, P.E.K.; investigation, P.E.K.;
resources, E.M.A. and J.W.C.; data curation, P.E.K.; writing—original draft,
P.E.K.; writing—review and editing, all; visualization, P.E.K.; supervision,
M.G., E.M.A., and J.W.C.; funding acquisition, E.M.A and J.W.C. All authors
have read and agreed to the published version of the manuscript.

Not applicable.

Not applicable. Not applicable.  Working on publicly available dataset.

###### Acknowledgements.

This work could not have been done without help from my friends in the
Michigan eXploration Lab (MXL) and Autonomous Aerospace Systems (A2sys) lab.
Specifically, thanks to Matthew Romano, Prashin Sharma, Akshay Mathur, Mark
Nail, Justin Schachter, Tran Anh Nguyen, and Paul Flanigen for their support
in fabricating my test vehicle, tuning the flight controller, assisting with
flight tests, developing/soldering PCBs for the magnetometers, and helping to
concisely communicating my work through the opening figure. Special thanks to
Mark Moldwin, Alex Hoffmann, and Brady Strabel from the University of Michigan
Magnetometer Laboratory (Mag Lab) for assistance with configuring the PNI
RM3100 magnetometer for this work. The authors declare no conflict of
interest. no

##  Appendix A Flight Test Name Reference

The, rather obtuse, naming convention for flight tests in this paper is meant
to serve the reader in associating our results with the relevant data in our
dataset. Recall that each flight test is denoted as “tY_XX” where Y refers to
the flight test series and XX is the two-digit ID of the flight test in that
series. The Table 5 associates each series numbers Y used in this paper with
their respective directories in our dataset.

Table 5: Data used in this paper.

Test | Dataset | Use In
---|---|---
Series (tY) | Directory | This Paper
t1 | 20220617_batteryAndMagTesting |  $S_{2}=2$cm
t2 | 20220622_batteryAndMagPartTwo |  $S_{2}=4$cm
t3 | 20220627_verifyingSufficientMagDist |  $S_{2}=6$cm
t4 | 20220707_8cm_magAndBatt |  $S_{2}=8$cm
t5 | 20220503_maggieAndEstibonFlights | ‘noisy’ versus ‘quiet’ PID controller gains
t6 | 20220901_researchData | Intermediate and compromise maps. Consistency metric.

References

## References

  * Kuevor et al. (2021) Kuevor, P.E.; Cutler, J.W.; Atkins, E.M.  Improving Attitude Estimation Using Gaussian-Process-Regression-Based Magnetic Field Maps.  Sensors 2021, 21, 6351.  https://doi.org/10.3390/s21196351.
  * Kok and Solin (2018) Kok, M.; Solin, A.  Scalable Magnetic Field SLAM in 3D Using Gaussian Process Maps.  2018 21st International Conference on Information Fusion (FUSION) 2018, pp. 1353–1360, [1804.01926].  https://doi.org/10.23919/ICIF.2018.8455789.
  * Vallivaara et al. (2011) Vallivaara, I.; Haverinen, J.; Kemppainen, A.; Roning, J.  Magnetic field-based SLAM method for solving the localization problem in mobile robot floor-cleaning task.  In Proceedings of the 2011 15th International Conference on Advanced Robotics (ICAR). IEEE, 2011, pp. 198–203.  https://doi.org/10.1109/ICAR.2011.6088632.
  * Vallivaara et al. (2010) Vallivaara, I.; Haverinen, J.; Kemppainen, A.; Roning, J.  Simultaneous localization and mapping using ambient magnetic field.  In Proceedings of the 2010 IEEE Conference on Multisensor Fusion and Integration. IEEE, 2010, pp. 14–19.  https://doi.org/10.1109/MFI.2010.5604465.
  * Akai and Ozaki (2015) Akai, N.; Ozaki, K.  Gaussian processes for magnetic map-based localization in large-scale indoor environments.  IEEE International Conference on Intelligent Robots and Systems 2015, 2015-December, 4459–4464.  https://doi.org/10.1109/IROS.2015.7354010.
  * Akai and Ozaki (2017) Akai, N.; Ozaki, K.  3D magnetic field mapping in large-scale indoor environment using measurement robot and Gaussian processes.  In Proceedings of the 2017 International Conference on Indoor Positioning and Indoor Navigation (IPIN). IEEE, 2017, Vol. 2017-Janua, pp. 1–7.  https://doi.org/10.1109/IPIN.2017.8115960.
  * Suksakulchai et al. (2000) Suksakulchai, S.; Thongchai, S.; Wilkes, D.; Kawamura, K.  Mobile robot localization using an electronic compass for corridor environment.  In Proceedings of the SMC 2000 Conference Proceedings. 2000 IEEE International Conference on Systems, Man and Cybernetics. ’Cybernetics Evolving to Systems, Humans, Organizations, and their Complex Interactions’ (Cat. No.00CH37166). IEEE, 2000, Vol. 5, pp. 3354–3359.  https://doi.org/10.1109/ICSMC.2000.886523.
  * Li et al. (2012) Li, B.; Gallagher, T.; Dempster, A.G.; Rizos, C.  How feasible is the use of magnetic field alone for indoor positioning?  In Proceedings of the 2012 International Conference on Indoor Positioning and Indoor Navigation (IPIN). IEEE, 2012, pp. 1–9.  https://doi.org/10.1109/IPIN.2012.6418880.
  * Haverinen and Kemppainen (2009) Haverinen, J.; Kemppainen, A.  Global indoor self-localization based on the ambient magnetic field.  Robotics and Autonomous Systems 2009, 57, 1028–1035.  https://doi.org/10.1016/j.robot.2009.07.018.
  * Wu et al. (2019) Wu, Z.; Wen, M.; Peng, G.; Tang, X.; Wang, D.  Magnetic-Assisted Initialization for Infrastructure-free Mobile Robot Localization.  In Proceedings of the 2019 IEEE International Conference on Cybernetics and Intelligent Systems (CIS) and IEEE Conference on Robotics, Automation and Mechatronics (RAM). IEEE, 2019, pp. 518–523, [1911.09313].  https://doi.org/10.1109/CIS-RAM47153.2019.9095809.
  * Robertson et al. (2013) Robertson, P.; Frassl, M.; Angermann, M.; Doniec, M.; Julian, B.J.; Puyol, M.G.; Khider, M.; Lichtenstern, M.; Bruno, L.  Simultaneous localization and mapping for pedestrians using distortions of the local magnetic field intensity in large indoor environments.  2013 International Conference on Indoor Positioning and Indoor Navigation, IPIN 2013 2013.  https://doi.org/10.1109/IPIN.2013.6817910.
  * Wang et al. (2017) Wang, X.; Zhang, C.; Liu, F.; Dong, Y.; Xu, X.  Exponentially weighted particle filter for simultaneous localization and mapping based on magnetic field measurements.  IEEE Transactions on Instrumentation and Measurement 2017, 66, 1658–1667.  https://doi.org/10.1109/TIM.2017.2664538.
  * Afzal et al. (2011) Afzal, M.H.; Renaudin, V.; Lachapelle, G.  Multi-magnetometer based perturbation mitigation for indoor orientation estimation.  Navigation, Journal of the Institute of Navigation 2011, 58, 279–292.  https://doi.org/10.1002/j.2161-4296.2011.tb02586.x.
  * Yin and Zhang (2018) Yin, G.; Zhang, L.  Magnetic heading compensation method based on magnetic interferential signal inversion.  Sensors and Actuators A: Physical 2018, 275, 1–10.  https://doi.org/10.1016/j.sna.2018.03.043.
  * Karimi et al. (2020) Karimi, M.; Babaians, E.; Oelsch, M.; Aykut, T.; Steinbach, E.  Skewed-redundant Hall-effect Magnetic Sensor Fusion for Perturbation-free Indoor Heading Estimation.  Proceedings - 4th IEEE International Conference on Robotic Computing, IRC 2020 2020, pp. 367–374.  https://doi.org/10.1109/IRC.2020.00064.
  * Hanley et al. (2021) Hanley, D.; De Oliveira, A.S.; Zhang, X.; Kim, D.H.; Wei, Y.; Bretl, T.  The Impact of Height on Indoor Positioning with Magnetic Fields.  IEEE Transactions on Instrumentation and Measurement 2021, 70.  https://doi.org/10.1109/TIM.2021.3059317.
  * Brzozowski et al. (2016) Brzozowski, B.; Kazmierczak, K.; Rochala, Z.; Wojda, M.; Wojtowicz, K.  A concept of UAV indoor navigation system based on magnetic field measurements.  3rd IEEE International Workshop on Metrology for Aerospace, MetroAeroSpace 2016 - Proceedings 2016, pp. 636–640.  https://doi.org/10.1109/MetroAeroSpace.2016.7573291.
  * Brzozowski and Kazmierczak (2017) Brzozowski, B.; Kazmierczak, K.  Magnetic field mapping as a support for UAV indoor navigation system.  In Proceedings of the 2017 IEEE International Workshop on Metrology for AeroSpace (MetroAeroSpace). IEEE, 2017, pp. 583–588.  https://doi.org/10.1109/MetroAeroSpace.2017.7999535.
  * Frassl et al. (2013) Frassl, M.; Angermann, M.; Lichtenstern, M.; Robertson, P.; Julian, B.J.; Doniec, M.  Magnetic maps of indoor environments for precise localization of legged and non-legged locomotion.  IEEE International Conference on Intelligent Robots and Systems 2013, pp. 913–920.  https://doi.org/10.1109/IROS.2013.6696459.
  * Liu et al. (2021) Liu, G.; Yu, B.; Huang, L.; Shi, L.; Gao, X.; He, L.  Human-Interactive Mapping Method for Indoor Magnetic Based on Low-Cost MARG Sensors.  IEEE Transactions on Instrumentation and Measurement 2021, 70.  https://doi.org/10.1109/TIM.2021.3052026.
  * Gozick et al. (2011) Gozick, B.; Subbu, K.P.; Dantu, R.; Maeshiro, T.  Magnetic maps for indoor navigation.  IEEE Transactions on Instrumentation and Measurement 2011, 60, 3883–3891.  https://doi.org/10.1109/TIM.2011.2147690.
  * Ayanoglu et al. (2018) Ayanoglu, A.; Schneider, D.M.; Eitel, B.  Crowdsourcing-Based Magnetic Map Generation for Indoor Localization.  In Proceedings of the 2018 International Conference on Indoor Positioning and Indoor Navigation (IPIN). IEEE, 2018, pp. 1–8.  https://doi.org/10.1109/IPIN.2018.8533832.
  * Almeida et al. (2021) Almeida, D.; Pedrosa, E.; Curado, F.  Magnetic Mapping for Robot Navigation in Indoor Environments.  2021 International Conference on Indoor Positioning and Indoor Navigation, IPIN 2021 2021.  https://doi.org/10.1109/IPIN51156.2021.9662528.
  * Zheng et al. (2021) Zheng, Y.; Li, S.; Xing, K.; Zhang, X.  Unmanned aerial vehicles for magnetic surveys: A review on platform selection and interference suppression.  Drones 2021, 5.  https://doi.org/10.3390/drones5030093.
  * Versteeg and McKay (2007) Versteeg, R.; McKay, M.  Feasibility Study for an Autonomous UAV -Magnetometer System – Final Report on SERDP SEED 1509:2206.  Technical Report September, Idaho National Laboratory, 2007.  https://doi.org/10.2172/923485.
  * Koyama et al. (2013) Koyama, T.; Kaneko, T.; Ohminato, T.; Yanagisawa, T.; Watanabe, A.; Takeo, M.  An aeromagnetic survey of Shinmoe-dake volcano, Kirishima, Japan, after the 2011 eruption using an unmanned autonomous helicopter.  Earth, Planets and Space 2013, 65, 657–666.  https://doi.org/10.5047/eps.2013.03.005.
  * Vasiljevic et al. (2022) Vasiljevic, G.; Martinovic, D.; Batos, M.; Bogdan, S.  Validation of two-wire power line UAV localization based on the magnetic field strength.  2022 International Conference on Unmanned Aircraft Systems, ICUAS 2022 2022, pp. 434–441.  https://doi.org/10.1109/ICUAS54217.2022.9836106.
  * Li et al. (2019) Li, Y.; Zahran, S.; Zhuang, Y.; Gao, Z.; Luo, Y.; He, Z.; Pei, L.; Chen, R.; El-Sheimy, N.  IMU/magnetometer/barometer/mass-flow sensor integrated indoor quadrotor UAV localization with robust velocity updates.  Remote Sensing 2019, 11, 1–22.  https://doi.org/10.3390/rs11070838.
  * Zahran et al. (2019) Zahran, S.; Moussa, A.M.; Sesay, A.B.; El-Sheimy, N.  A New Velocity Meter Based on Hall Effect Sensors for UAV Indoor Navigation.  IEEE Sensors Journal 2019, 19, 3067–3076.  https://doi.org/10.1109/JSEN.2018.2890094.
  * Tiemann et al. (2018) Tiemann, J.; Ramsey, A.; Wietfeld, C.  Enhanced UAV indoor navigation through SLAM-Augmented UWB Localization.  2018 IEEE International Conference on Communications Workshops, ICC Workshops 2018 - Proceedings 2018, pp. 1–6.  https://doi.org/10.1109/ICCW.2018.8403539.
  * Niu et al. (2021) Niu, G.; Zhang, J.; Guo, S.; Pun, M.O.; Chen, C.S.  UAV-Enabled 3D Indoor Positioning and Navigation Based on VLC.  In Proceedings of the ICC 2021 - IEEE International Conference on Communications. IEEE, 2021, pp. 1–6.  https://doi.org/10.1109/ICC42927.2021.9500633.
  * Lipovský et al. (2021) Lipovský, P.; Draganová, K.; Novotňák, J.; Szőke, Z.; Fiľko, M.  Indoor mapping of magnetic fields using uav equipped with fluxgate magnetometer.  Sensors 2021, 21.  https://doi.org/10.3390/s21124191.
  * Blažek et al. (2018) Blažek, J.; Lipovský, P.; Heško, F.; Repčík, D.  Electromagnetic image of small UAV in very low frequency range.  Journal of Electrical Engineering 2018, 69, 438–441.  https://doi.org/10.2478/jee-2018-0069.
  * Rasmussen and Williams (2006) Rasmussen, C.E.; Williams, C.K.I.  Gaussian Processes for Machine Learning; MIT Press, 2006.
  * Springmann and Cutler (2012) Springmann, J.C.; Cutler, J.W.  Attitude-Independent Magnetometer Calibration with Time-Varying Bias.  Journal of Guidance, Control, and Dynamics 2012, 35, 1080–1088.  https://doi.org/10.2514/1.56726.
  * Wu et al. (2020) Wu, H.; Pei, X.; Li, J.; Gao, H.; Bai, Y.  An improved magnetometer calibration and compensation method based on Levenberg–Marquardt algorithm for multi-rotor unmanned aerial vehicle.  Measurement and Control 2020, 53, 276–286.  https://doi.org/10.1177/0020294019890627.
  * Team and Survey (2020) Team, N.G.M.; Survey, B.G.  World Magnetic Model 2020, 2020.  https://doi.org/10.25921/11v3-da71.
  * Regoli et al. (2018) Regoli, L.H.; Moldwin, M.B.; Pellioni, M.; Bronner, B.; Hite, K.; Sheinker, A.; Ponder, B.M.  Investigation of a low-cost magneto-inductive magnetometer for space science applications.  Geoscientific Instrumentation, Methods and Data Systems 2018, 7, 129–142.  https://doi.org/10.5194/gi-7-129-2018.
  * Romano et al. (2022) Romano, M.; Ye, A.; Pye, J.; Atkins, E.  Cooperative Multilift Slung Load Transportation Using Haptic Admittance Control Guidance.  Journal of Guidance, Control, and Dynamics 2022, pp. 1–14.  https://doi.org/10.2514/1.G006587.
  * Le Grand and Thrun (2012) Le Grand, E.; Thrun, S.  3-Axis magnetic field mapping and fusion for indoor localization.  In Proceedings of the 2012 IEEE International Conference on Multisensor Fusion and Integration for Intelligent Systems (MFI). IEEE, 2012, pp. 358–364.  https://doi.org/10.1109/MFI.2012.6343024.
  * Vandermeulen et al. (2013) Vandermeulen, D.; Vercauteren, C.; Weyn, M.  Indoor localization Using a Magnetic Flux Density Map of a Building Feasibility study of geomagnetic indoor localization.  AMBIENT 2013 : The Third International Conference on Ambient Computing, Applications, Services and Technologies Magnetic 2013, pp. 42–49.
  * Solin et al. (2015) Solin, A.; Kok, M.; Wahlström, N.; Schön, T.B.; Särkkä, S.  Modeling and interpolation of the ambient magnetic field by Gaussian processes.  IEEE Transactions on Robotics 2015, 34, 1112–1127, [1509.04634].  https://doi.org/10.1109/TRO.2018.2830326.
  * Michel et al. (2018) Michel, T.; Genevès, P.; Fourati, H.; Layaïda, N.  Attitude estimation for indoor navigation and augmented reality with smartphones.  Pervasive and Mobile Computing 2018, 46, 96–121.  https://doi.org/10.1016/j.pmcj.2018.03.004.