# Efficient Extrinsic Self-Calibration of Multiple IMUs using Measurement
Subset Selection

Jongwon Lee1, David Hanley2, and Timothy Bretl1 1Jongwon Lee and Timothy Bretl
are with the Department of Aerospace Engineering, University of Illinois
Urbana-Champaign, Urbana, IL 61801, USA (Email: {jongwon5,
tbretl}@illinois.edu).2David Hanley is with the School of Informatics and the
Translational Healthcare Technologies Group in Centre for Inflammation
Research, Institute for Regeneration and Repair, University of Edinburgh, UK
(Email: dhanley@ed.ac.uk).This work was supported by the NASA Grant No.
STTR-80NSSC20C0020.

###### Abstract

This paper addresses the problem of choosing a sparse subset of measurements
for quick calibration parameter estimation. A standard solution to this is
selecting a measurement only if its utility—the difference between posterior
(with the measurement) and prior information (without the measurement)—exceeds
some threshold. Theoretically, utility, a function of the parameter estimate,
should be evaluated at the estimate obtained with all measurements selected so
far, hence necessitating a recalibration with each new measurement. However,
we hypothesize that utility is insensitive to changes in the parameter
estimate for many systems of interest, suggesting that evaluating utility at
some initial parameter guess would yield equivalent results in practice. We
provide evidence supporting this hypothesis for extrinsic calibration of
multiple inertial measurement units (IMUs), showing the reduction in
calibration time by two orders of magnitude by forgoing recalibration for each
measurement.

##  I Introduction

Inertial measurement units (IMUs) are a common component in mobile robot
navigation, comprising accelerometers and gyroscopes to measure specific
forces and angular rates along three orthogonal axes. While single IMUs are
commonly used, deploying multiple IMUs can enhance measurement accuracy,
bandwidth, and fault tolerance without increasing the overall size, weight,
power, and cost [1, 2, 3, 4, 5, 6].

To realize these benefits, accurate extrinsic calibration is crucial to
determine the relative pose of each IMU to others. Various extrinsic
calibration methods for multi-IMU systems have been developed, including those
utilizing prescribed trajectories [7, 8, 9] and aiding sensors like cameras
[10]. Notably, self-calibration methods [11, 12, 13], which rely solely on IMU
measurements themselves, offer significant benefits, especially in scenarios
requiring recalibration due to intentional or unintentional changes in sensor
configuration during robot operation, such as part loosening or thermal
expansion.

Self-calibration, not bound to a predetermined trajectory, naturally
experiences measurements with varying “utility”. In cases like calibrating
multi-IMU systems on non-rotating bodies such as spacecraft in orbit or cars
on straight paths, much of the data can be low in utility. This is especially
critical for optimization-based calibration, favored for its accuracy compared
to filter-based methods, as it involves processing a significant amount of
low-utility (uninformative) data alongside with high-utility (informative)
data during optimization. Therefore, identifying subsets of measurements with
high utility would offer benefits such as reduced runtime and memory usage for
optimization.

In this study, we present an approach to multi-IMU extrinsic calibration by
selecting informative measurement subsets with high utility, significantly
streamlining our previous work [13]. This selection is primarily based on
method from existing studies [14, 15] but is executed more efficiently.
Specifically, we hypothesize that the measure of utility is insensitive to
changes in the parameter estimate and claim that equivalent calibration
results are obtained by evaluating utility only at some initial guess of the
calibration parameters. A valuable consequence of this is eliminating the need
for recalibration with each new measurement, thereby reducing total
computation time for calibration.

Our paper is structured as follows: Section II reviews existing literature on
informative measurement selection for self-calibration. Section III describes
an approach with our hypothesis for efficient measurement subset selection.
Section IV delineates the solution approach for multi-IMU extrinsic self-
calibration using these subsets, building on our prior work [13]. We then
proceed to assess its runtime and accuracy against benchmarks in simulation
(Section V) and using real-world data (Section VI). Section VII summarizes our
findings, provides conclusions, and discuss the implications of this work.

##  II Related Works

Selecting measurements with high utility for self-calibration involves
optimizing the Fisher information matrix $\mathbf{\bm{\mathcal{I}}}$, whose
inverse sets the Cramer-Rao lower bound [16] for parameter estimation. Imagine
a full measurement set $\mathcal{D}$ divided into $L$ segments
${\mathcal{D}^{1},\dotsc,\mathcal{D}^{L}}$, with each segment
$\mathcal{D}^{l}$ containing measurements across $K$ consecutive timesteps.
The goal is to identify a subset of $\mathcal{D}$

$\displaystyle\underset{\mathcal{D}^{info}\subset\mathcal{D},\,\text{card}(\mathcal{D}^{info})=\kappa}{\operatorname*{arg\,max}}$
$\displaystyle f\left[\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info})\right]$
(original) subject to $\displaystyle\kappa<\text{card}(\mathcal{D}),$

where $f$ maps the input matrix to a scalar value (termed “information” in
what follows), with $\text{card}(\cdot)$ denoting set cardinality and $\kappa$
a predetermined subset size. This combinatorial optimization problem, aimed at
maximizing information, is typically addressed through two alternative ways in
the existing literature. One approach involves choosing segments with the
highest information:

  1. 1.

Evaluate $f[\mathbf{\bm{\mathcal{I}}}(\cdot)]$ for segments $\mathcal{D}^{1}$
through $\mathcal{D}^{L}$.

  2. 2.

Choose $M$ segments showing the largest $f[\mathbf{\bm{\mathcal{I}}}(\cdot)]$.

(M-largest selection)

Schneider et al. [17] and Lv et al. [18] used this approach in camera-IMU and
LiDAR-IMU calibration, respectively. The second approach uses a greedy
algorithm, which selectively adds candidate segments $\mathcal{D}^{new}$ in
the informative subset $\mathcal{D}^{info}$ if their utility exceeds threshold
$\lambda$:

  1. 1.

Initialize $\mathcal{D}^{info}$ as an empty set.

  2. 2.

Iterate over candidate segments $\mathcal{D}^{new}$ from $\mathcal{D}^{1}$ to
$\mathcal{D}^{L}$.

  3. 3.

Add $\mathcal{D}^{new}$ in $\mathcal{D}^{info}$ if
$f[\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})]-f[\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info})]>\lambda$.

  4. 4.

Continue until all segments are processed.

(greedy algorithm)

Maye et al. initially applied this approach in rangefinder-odometer [14] and
later in camera-IMU calibration [15].

Both approaches succeed in calibration only if they correctly identify
informative subsets. We show in Section V-B an edge case where M-largest
selection fails, prompting our preference for the greedy algorithm.

However, the usefulness of the greedy algorithm depends on the selection
efficiency. We show that the original greedy algorithm by Maye et al. [14, 15]
may not significantly improve or might even worsen runtime for our multi-IMU
calibration scenario in concern. This led us to seek improvements detailed in
Section III, where we apply our hypothesis for a more efficient execution of
the greedy algorithm.

##  III Efficient Identification of Informative Measurement Subsets for Self-
Calibration

Algorithm 1 Pseudocode for calibration via the original greedy algorithm [14,
15]

1:Input:  uncalibrated parameters $\mathbf{\bm{\theta}}^{0}$,

2:Output:  calibrated parameters $\mathbf{\bm{\theta}}^{\ast}$

3:$\mathcal{D}^{info}\leftarrow\varnothing$

4:$\widehat{\mathbf{\bm{\theta}}}^{-}\leftarrow\mathbf{\bm{\theta}}^{0}$

5:for $new\leftarrow 1$ to $L$ do

6:
$\widehat{\mathbf{\bm{\theta}}}^{+}\leftarrow\texttt{Calibrate}(\mathcal{D}^{info},\mathcal{D}^{new};\widehat{\mathbf{\bm{\theta}}}^{-})$

7:
$\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})\leftarrow\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\widehat{\mathbf{\bm{\theta}}}^{+}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\widehat{\mathbf{\bm{\theta}}}^{+}}+\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\widehat{\mathbf{\bm{\theta}}}^{+}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\widehat{\mathbf{\bm{\theta}}}^{+}}$

8:
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}\leftarrow\big{[}\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})^{-1}\big{]}_{\mathbf{\bm{\Theta}}\mathbf{\bm{\Theta}}}$

9: if
$\frac{1}{2}\textrm{log}\frac{\big{|}\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info}}\big{|}}{\big{|}\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}\big{|}}>\lambda$
or $\mathcal{D}^{info}==\varnothing$ then

10: $\mathcal{D}^{info}\leftarrow\mathcal{D}^{info}\cup\mathcal{D}^{new}$

11:
$|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info}}|\leftarrow|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}|$

12:
$\widehat{\mathbf{\bm{\theta}}}^{-}\leftarrow\widehat{\mathbf{\bm{\theta}}}^{+}$

13: end if

14:end for

15:$\mathbf{\bm{\theta}}^{\ast}\leftarrow\widehat{\mathbf{\bm{\theta}}}^{-}$

16:return $\mathbf{\bm{\theta}}^{\ast}$

Algorithm 2 Pseudocode for calibration via greedy algorithm with utility
evaluation at initial calibration parameters

1:Input:  uncalibrated parameters $\mathbf{\bm{\theta}}^{0}$,

2:Output:  calibrated parameters $\mathbf{\bm{\theta}}^{\ast}$

3:$\mathcal{D}^{info}\leftarrow\varnothing$

4:for $new\leftarrow 1$ to $L$ do

5:
$\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})\leftarrow\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}+\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\theta}}^{0}}$

6:
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}\leftarrow\big{[}\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})^{-1}\big{]}_{\mathbf{\bm{\Theta}}\mathbf{\bm{\Theta}}}$

7: if
$\frac{1}{2}\textrm{log}\frac{\big{|}\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info}}\big{|}}{\big{|}\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}\big{|}}>\lambda$
or $\mathcal{D}^{info}==\varnothing$ then

8: $\mathcal{D}^{info}\leftarrow\mathcal{D}^{info}\cup\mathcal{D}^{new}$ (Keep
$\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}$
as well)

9:
$|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info}}|\leftarrow|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}|$

10: end if

11:end for

12:$\mathbf{\bm{\theta}}^{\ast}\leftarrow\texttt{Calibrate}(\mathcal{D}^{info};\mathbf{\bm{\theta}}^{0})$

13:return $\mathbf{\bm{\theta}}^{\ast}$

In this section, we first present the greedy algorithm by Maye et al. [14, 15]
(Alg. 1), then present its adaptation, by evaluating utility at initial
calibration parameters (Alg. 2). We provide a theoretical analysis showing its
computational complexity reduction, expressed in Big O notation.

###  III-A Original Greedy Algorithm

We provide details of Alg. 1 here. Upon receiving new measurements
$\mathcal{D}^{new}$, the algorithm obtains the calibration parameter estimate
$\mathbf{\bm{\widehat{\theta}}}^{+}$ by minimizing the squared error sum
associated with both $\mathcal{D}^{info}$ and $\mathcal{D}^{new}$, starting
from an initial estimate $\mathbf{\bm{\widehat{\theta}}}^{-}$, through the
optimization process (Line 6, Alg. 1). It then evaluates the Jacobian for both
sets at the current estimate $\mathbf{\bm{\widehat{\theta}}}^{+}$
($\mathbf{\bm{J}}(\mathcal{D}^{info},\mathcal{D}^{new})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}=\begin{bmatrix}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}^{T}&\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}^{T}\end{bmatrix}^{T}$)
and computes the Fisher information matrix:
$\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})=\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}+\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\widehat{\theta}}}^{+}}$
(Line 7, Alg. 1). Next, the block entries corresponding to
$\mathbf{\bm{\Theta}}$ from its inverse are taken:
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}=\big{[}\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})^{-1}\big{]}_{\mathbf{\bm{\Theta}}\mathbf{\bm{\Theta}}}$
(Line 8, Alg. 1). If the utility of $\mathcal{D}^{new}$ — the difference
between
$-1/2\cdot\log|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}|$
and
$-1/2\cdot\log|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info}}|$
(obtained from the previous iteration) — exceeds $\lambda$, then
$\mathcal{D}^{new}$ is considered a high-utility segment (Line 9, Alg. 1),
prompting the parameter updates (Lines 10-12, Alg. 1). Note that this process
iteratively refines calibration.

###  III-B Greedy Algorithm with Utility Evaluation at Initial Calibration
Parameters

Alg. 2 distinguishes itself from Alg. 1 by eliminating parameter recalibration
in its loop. This is based on the hypothesis that utility evaluation for
measurements is unaffected by the specific choice of $\mathbf{\bm{\theta}}$,
thus not significantly altering the choice of $\mathcal{D}^{info}$ or the
calibration results.

Upon receiving new measurements $\mathcal{D}^{new}$, Alg. 2 evaluates the
Jacobian _only_ for $\mathcal{D}^{new}$ using uncalibrated parameters
$\mathbf{\bm{\theta}}^{0}$, computing the Fisher information matrix:
$\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})=\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}+\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\theta}}^{0}}$
(Line 5, Alg. 2). The repetitive calculation of
$\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}$
is omitted in each iteration, by using the result from the preceding iteration
(Line 8, Alg. 2). Then, the block entries corresponding to
$\mathbf{\bm{\Theta}}$ from its inverse are taken:
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}=\big{[}\mathbf{\bm{\mathcal{I}}}(\mathcal{D}^{info},\mathcal{D}^{new})^{-1}\big{]}_{\mathbf{\bm{\Theta}}\mathbf{\bm{\Theta}}}$
(Line 6, Alg. 2). If the utility of $\mathcal{D}^{new}$ — the difference
between
$-1/2\cdot\log|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info},\mathcal{D}^{new}}|$
and
$-1/2\cdot\log|\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\Theta}}|\mathcal{D}^{info}}|$
(obtained from the previous iteration) — exceeds $\lambda$, then
$\mathcal{D}^{new}$ is considered a high-utility segment (Line 7, Alg. 2),
prompting the parameter updates and retention of
$\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}^{T}\mathbf{\bm{J}}(\mathcal{D}^{info})|_{\mathbf{\bm{\theta}}^{0}}$
for future iterations (Lines 8-9, Alg. 2). Finally, with $\mathcal{D}^{info}$
identified, the calibration is conducted by solving a nonlinear least-squares
problem (Line 12, Alg. 2).

###  III-C Comparison of Time Complexity

Consequently, Alg. 2 yields reduced time complexity compared to Alg. 1, as
shown in Table II. Detailed derivations are provided in Appendices A and B,
respectively.

Alg. 2 is notably effective when the high-utility subset selection
$\mathcal{D}^{info}$ is not highly sensitive to variations in
$\mathbf{\bm{\theta}}$. To validate this, we conduct a sensitivity analysis on
the impact of $\mathbf{\bm{\theta}}$ deviations in measured information both
in simulations (Section V-C) and hardware experiments (Section VI-B). We also
show that calibration results between Alg. 1 and Alg. 2 are comparably
insignificant in both the simulation (Section V) and hardware experiment
sections (Section VI).

TABLE I: Big O Comparison: Original Greedy vs. Greedy with Utility Evaluation
at Initial Calibration Parameters

| Original | Init-Param
---|---|---
Best | $\mathcal{O}(LKp^{2})$ | $\mathcal{O}(LKp^{2})$
Worst | $\mathcal{O}(L^{2}Kp^{2})$ | $\mathcal{O}(LKp^{2})$
|  |

TABLE II: *

* $L$: the number of total segments; $K$: the number of timesteps (i.e., data) in each segment; $p$: the number of parameters being estimated

##  IV Self-Calibration for Multiple IMUs

This section presents the solution approach for multi-IMU extrinsic self-
calibration with sparse measurement subsets, beginning with a problem overview
from Section IV-A to Section IV-C, which serves as a review of our previous
work [13]. We detail a method for estimating each IMU’s relative orientation
in Section IV-D, crucial for Alg. 2’s validity, as shown in Section V-C and
VI-B. Section IV-E discusses residuals for evaluating the Fisher information
matrix and formulating the calibration’s nonlinear least-squares problem.
Section IV-F introduces additional residuals for constraining the relationship
between the sparse measurement subsets, followed by the calibration solution
approach.

###  IV-A Notation

We denote the rotation matrix that represents frame $\mathcal{F}_{A}$’s
orientation in frame $\mathcal{F}_{B}$’s coordinates by
${}^{B}_{A}\mathbf{\bm{R}}\>\in\>SO(3)$. The operator $\mathbf{\bm{C}}(\cdot)$
converts a unit quaternion ${}^{B}_{A}\mathbf{\bm{q}}$ to its corresponding
rotation matrix ${}^{B}_{A}\mathbf{\bm{R}}$. Coordinate transformations are
thus defined:

${}^{B}\mathbf{\bm{p}}={}^{B}_{A}\mathbf{\bm{R}}{}^{A}\mathbf{\bm{p}}=\mathbf{\bm{C}}\big{(}{}^{B}_{A}\mathbf{\bm{q}}\big{)}{}^{A}\mathbf{\bm{p}}.$

We distinguish quantities using these symbols:

  * •

a tilde $\widetilde{(\,\cdot\,)}$ denotes a sensor measurement,

  * •

a hat $\widehat{(\,\cdot\,)}$ denotes an estimate,

  * •

a zero superscript $(\,\cdot\,)^{0}$ denotes an initial guess (pre-
optimization),

  * •

and an asterisk superscript $(\,\cdot\,)^{*}$ denotes a final value (post-
optimization).

###  IV-B Inertial Sensor Model

We model accelerometer measurements from each IMU as

$^{I}\widetilde{\mathbf{\bm{a}}}={}^{I}\mathbf{\bm{a}}_{WI}-{}^{I}\mathbf{\bm{g}}+\mathbf{\bm{b}}_{a}+\mathbf{\bm{n}}_{a}$
(1)

where

  * •

${}^{I}\mathbf{\bm{a}}_{WI}$ is the IMU’s linear acceleration relative to the
world frame, expressed in the IMU frame,

  * •

${}^{I}\mathbf{\bm{g}}$ is gravitational acceleration in the IMU frame,

  * •

$\mathbf{\bm{b}}_{a}$ is a time-varying bias modeled by a random walk

$\mathbf{\bm{b}}_{a,k+1}-\mathbf{\bm{b}}_{a,k}\sim\sigma_{\mathbf{\bm{b}}_{a}}\sqrt{\Delta
t}\cdot\mathcal{N}\left(\mathbf{\bm{0}},\,\mathbf{\bm{1}}\right),$

  * •

and $\mathbf{\bm{n}}_{a}$ is stochastic noise

$\mathbf{\bm{n}}_{a}\sim\sigma_{a}/\sqrt{\Delta
t}\cdot\mathcal{N}\left(\mathbf{\bm{0}},\,\mathbf{\bm{1}}\right)$

with $\Delta t$ as the IMU’s sampling interval. We model gyroscope
measurements from each IMU as

$^{g}\widetilde{\mathbf{\bm{\omega}}}=\mathbf{\bm{C}}\big{(}{}^{g}_{I}\mathbf{\bm{q}}\big{)}{}^{I}\mathbf{\bm{\omega}}_{WI}+\mathbf{\bm{b}}_{g}+\mathbf{\bm{n}}_{g}$
(2)

where

  * •

$\mathbf{\bm{C}}(^{g}_{I}\mathbf{\bm{q}})$ is the rotation matrix (written in
terms of the corresponding quaternion) that models the IMU frame’s orientation
relative to the gyroscope frame, accounting for gyroscope misalignment,

  * •

${}^{I}\mathbf{\bm{\omega}}_{WI}$ is the IMU’s angular velocity relative to
the world frame, expressed in the IMU frame,

  * •

$\mathbf{\bm{b}}_{g}$ is a time-varying bias modeled a random walk

$\mathbf{\bm{b}}_{g,k+1}-\mathbf{\bm{b}}_{g,k}\sim\sigma_{\mathbf{\bm{b}}_{g}}\sqrt{\Delta
t}\cdot\mathcal{N}\left(\mathbf{\bm{0}},\,\mathbf{\bm{1}}\right),$

  * •

and $\mathbf{\bm{n}}_{g}$ is stochastic noise

$\mathbf{\bm{n}}_{g}\sim\sigma_{g}/\sqrt{\Delta
t}\cdot\mathcal{N}\left(\mathbf{\bm{0}},\,\mathbf{\bm{1}}\right).$

We assume that noise densities $\sigma_{a}$, $\sigma_{g}$ and bias
instabilities $\sigma_{\mathbf{\bm{b}}_{a}}$, $\sigma_{\mathbf{\bm{b}}_{g}}$
in continuous-time are identified for all sensors before extrinsic
calibration. Additionally, we assume intrinsic parameters like scale factors
and axis non-orthogonality are calibrated, and all IMUs are temporally
synchronized.

###  IV-C Problem Statement (Multi-IMU Extrinsic Calibration)

Supposing $N+1$ IMUs whose frames are indexed as $I_{0},\,\dotsc,\,I_{N}$, we
define:

  * •

${}^{I_{n}}\mathbf{\bm{p}}_{I_{n}I_{0}}:=\mathbf{\bm{p}}_{{I}_{n}}$ for the
frame $I_{0}$’s position relative to frame $I_{n}$, expressed in frame
$I_{n}$,

  * •

${}^{I_{n}}_{I_{0}}\mathbf{\bm{q}}:=\mathbf{\bm{q}}_{I_{n}}$ for the
quaternion describing frame $I_{0}$’s orientation relative to frame $I_{n}$,
assuming IMU frames are accelerometer-aligned,

  * •

and ${}^{I_{n}}_{g_{n}}\mathbf{\bm{q}}:=\mathbf{\bm{q}}_{g_{n}}$ for the
quaternion describing frame $g_{n}$’s orientation relative to frame $I_{n}$
(gyroscope misalignment).

Our calibration aims to estimate extrinsic parameters
$\mathbf{\bm{p}}_{{I}_{n}}$ and $\mathbf{\bm{q}}_{{I}_{n}}$ for each IMU
$n\in\\{1,\dotsc,N\\}$ and $\mathbf{\bm{q}}_{g_{n}}$ for each IMU
$n\in\\{0,\dotsc,N\\}$, given measurements
${}^{I_{n}}\widetilde{\mathbf{\bm{a}}}_{k},\,^{g_{n}}\widetilde{\mathbf{\bm{\omega}}}_{k}$
at each time $k$ for each IMU $n\in\\{0,\dotsc,N\\}$. Additionally, estimating
time-varying biases $\mathbf{\bm{b}}_{a_{n},k},\,\mathbf{\bm{b}}_{g_{n},k}$ at
each time $k$ for each IMU $n\in\\{0,\dotsc,N\\}$, as well as the angular
acceleration of the base IMU
${}^{I_{0}}\mathbf{\bm{\alpha}}_{I_{0},k}:=^{I_{0}}\mathbf{\bm{\alpha}}_{k}$
at each time $k$, proves useful in refining sensor parameter estimates in our
previous study [13].

###  IV-D Relative Orientation Initialization

Prior to subset selection in Alg. 2, we initialize each IMU’s relative
orientation to the base IMU, ${}^{I_{0}}_{I_{n}}\mathbf{\bm{q}}$, using Yang
et al.’s method [19]. This preliminary step enhances Alg. 2’s validity by
providing plausible IMU orientation guesses.

Assuming angular velocity $\mathbf{\bm{\omega}}$ remains constant between
$t_{k}$ and $t_{k+1}$, with $\Delta t=t_{k+1}-t_{k}$, the orientation change
in the world frame expressed in the IMU frame is

${}^{I}_{W}\mathbf{\bm{q}}(t_{k+1})=\begin{bmatrix}\frac{\mathbf{\bm{\omega}}}{|\mathbf{\bm{\omega}}|}\sin\left(\frac{|\mathbf{\bm{\omega}}|}{2}\Delta
t\right)\\\ \cos\left(\frac{|\mathbf{\bm{\omega}}|}{2}\Delta
t\right)\end{bmatrix}\;\otimes\;{}^{I}_{W}\mathbf{\bm{q}}(t_{k}).$

Hence, the rotation from $\mathcal{F}_{I(t_{k})}$ to
$\mathcal{F}_{I(t_{k+1})}$ is

${}^{I(t_{k+1})}_{I(t_{k})}\mathbf{\bm{q}}=\begin{bmatrix}\frac{\mathbf{\bm{\omega}}}{|\mathbf{\bm{\omega}}|}\sin\left(\frac{|\mathbf{\bm{\omega}}|}{2}\Delta
t\right)\\\ \cos\left(\frac{|\mathbf{\bm{\omega}}|}{2}\Delta
t\right)\end{bmatrix}.$

Supposing that we have ${}^{I_{0}(t_{k+1})}_{I_{0}(t_{k})}\mathbf{\bm{q}}$ and
${}^{I_{n}(t_{k+1})}_{I_{n}(t_{k})}\mathbf{\bm{q}}$ from
${}^{I_{0}}\widetilde{\mathbf{\bm{\omega}}}_{k}$ and
${}^{I_{n}}\widetilde{\mathbf{\bm{\omega}}}_{k}$ respectively, the following
relationship holds:

${}^{I_{0}(t_{k+1})}_{I_{0}(t_{k})}\mathbf{\bm{q}}\,\otimes\,^{I_{0}}_{I_{n}}\mathbf{\bm{q}}\;=\;^{I_{0}}_{I_{n}}\mathbf{\bm{q}}\,\otimes\,^{I_{n}(t_{k+1})}_{I_{n}(t_{k})}\mathbf{\bm{q}}.$

Or equivalently,

$\mathcal{L}(^{I_{0}(t_{k+1})}_{I_{0}(t_{k})}\mathbf{\bm{q}})\,\cdot\,^{I_{0}}_{I_{n}}\mathbf{\bm{q}}=\mathcal{R}(^{I_{n}(t_{k+1})}_{I_{n}(t_{k})}\mathbf{\bm{q}})\,\cdot\,^{I_{0}}_{I_{n}}\mathbf{\bm{q}}$

where

$\displaystyle\mathcal{L}(\mathbf{\bm{q}})$
$\displaystyle=\begin{bmatrix}q_{w}\mathbf{\bm{I}}_{3}+\lfloor\mathbf{\bm{q}}_{xyz}\rfloor&\mathbf{\bm{q}}_{xyz}\\\
-\mathbf{\bm{q}}_{xyz}^{T}&q_{w}\end{bmatrix},$
$\displaystyle\mathcal{R}(\mathbf{\bm{q}})$
$\displaystyle=\begin{bmatrix}q_{w}\mathbf{\bm{I}}_{3}-\lfloor\mathbf{\bm{q}}_{xyz}\rfloor&\mathbf{\bm{q}}_{xyz}\\\
-\mathbf{\bm{q}}_{xyz}^{T}&q_{w}\end{bmatrix}.$

Rearranging the equation yields

$\left(\mathcal{L}(^{I_{0}(t_{k+1})}_{I_{0}(t_{k})}\mathbf{\bm{q}})-\mathcal{R}(^{I_{n}(t_{k+1})}_{I_{n}(t_{k})}\mathbf{\bm{q}})\right)\,\cdot\,^{I_{0}}_{I_{n}}\mathbf{\bm{q}}=\mathbf{\bm{0}}_{4\times
1}.$

Stacking the difference of 4-by-4 matrices $\mathcal{L}-\mathcal{R}$ for
various timesteps $k$ creates a large matrix
$\mathbf{\bm{A}}\in\mathbb{R}^{4k\times 4}$, leading to
$\mathbf{\bm{A}}\,\cdot\,^{I_{0}}_{I_{n}}\mathbf{\bm{q}}=\mathbf{\bm{0}}_{4\times
1}$. We compute ${}^{I_{0}}_{I_{n}}\mathbf{\bm{q}}$ by finding the right-
singular vector of $\mathbf{\bm{A}}$ with the smallest singular value.

###  IV-E Intra-Segment Residuals

For each segment $\mathcal{D}^{new}$ with $N+1$ IMU measurements over $K$
timesteps, we define four types of residuals, $\mathbf{\bm{r}}_{a}$,
$\mathbf{\bm{r}}_{g}$, $\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}$, and
$\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}$. Residual $\mathbf{\bm{r}}_{a}$
relates the $n$th accelerometer’s measurement
(${}^{I_{n}}\widetilde{\mathbf{\bm{a}}}_{k}$) to that of base accelerometer’s
(${}^{I_{0}}\widetilde{\mathbf{\bm{a}}}_{k}$) corrected for bias
$\mathbf{\bm{b}}_{a_{0},k}$:

$\mathbf{\bm{r}}_{a}={}^{I_{0}}\widehat{\mathbf{\bm{a}}}_{k}-({}^{I_{0}}\widetilde{\mathbf{\bm{a}}}_{k}-\mathbf{\bm{b}}_{a_{0},k}).$

${}^{I_{0}}\widehat{\mathbf{\bm{a}}}_{k}$ is derived by transforming
${}^{I_{n}}\widetilde{\mathbf{\bm{a}}}_{k}$ to the base IMU’s frame:

$\displaystyle{}^{I_{0}}\widehat{\mathbf{\bm{a}}}_{k}=\;$
$\displaystyle{}^{I_{n}}_{I_{0}}\mathbf{\bm{R}}^{-1}\Big{\\{}({}^{I_{n}}\widetilde{\mathbf{\bm{a}}}_{k}-\mathbf{\bm{b}}_{a_{n},k})+$
$\displaystyle\lfloor{}^{I_{n}}_{g_{n}}\mathbf{\bm{R}}(^{g_{n}}\widetilde{\mathbf{\bm{\omega}}}_{k}-\mathbf{\bm{b}}_{g_{n},k})\rfloor^{2}\,\mathbf{\bm{p}}_{I_{n}}+\lfloor{}^{I_{n}}_{I_{0}}\mathbf{\bm{R}}^{I_{0}}{\mathbf{\bm{\alpha}}}_{k}\rfloor\,\mathbf{\bm{p}}_{I_{n}}\Big{\\}},$

where
${}^{I_{n}}_{I_{0}}\mathbf{\bm{R}}=\mathbf{\bm{C}}\big{(}\mathbf{\bm{q}}_{I_{n}}\big{)}$,
${}^{I_{n}}_{g_{n}}\mathbf{\bm{R}}=\mathbf{\bm{C}}\big{(}\mathbf{\bm{q}}_{g_{n}}\big{)}$,
and $\lfloor\cdot\rfloor$ mapping a cross product into skew-symmetric matrix
form. Residual $\mathbf{\bm{r}}_{g}$ relates the $n$th gyroscope’s measurement
(${}^{I_{n}}\widetilde{\mathbf{\bm{\omega}}}_{k}$) to that of base gyroscope’s
(${}^{I_{0}}\widetilde{\mathbf{\bm{\omega}}}_{k}$) corrected for biases
$\mathbf{\bm{b}}_{g_{n},k}$, $\mathbf{\bm{b}}_{g_{0},k}$:

$\mathbf{\bm{r}}_{g}={}^{I_{n}}_{I_{0}}\mathbf{\bm{R}}^{-1}{}^{I_{n}}_{g_{n}}\mathbf{\bm{R}}({}^{g_{n}}\widetilde{\mathbf{\bm{\omega}}}_{k}-\mathbf{\bm{b}}_{g_{n},k})-{}^{I_{0}}_{g_{0}}\mathbf{\bm{R}}({}^{g_{0}}\widetilde{\mathbf{\bm{\omega}}}_{k}-\mathbf{\bm{b}}_{g_{0},k}),$

where
${}^{I_{n}}_{I_{0}}\mathbf{\bm{R}}=\mathbf{\bm{C}}\big{(}\mathbf{\bm{q}}_{I_{n}}\big{)}$,
${}^{I_{n}}_{g_{n}}\mathbf{\bm{R}}=\mathbf{\bm{C}}\big{(}\mathbf{\bm{q}}_{g_{n}}\big{)}$,
and
${}^{I_{0}}_{g_{0}}\mathbf{\bm{R}}=\mathbf{\bm{C}}\big{(}\mathbf{\bm{q}}_{g_{0}}\big{)}$.
This relation holds as rotational motion is uniform across all points of a
rigid body. Sensor measurements corrected for bias introduce uncertainty due
to zero-mean Gaussian white noise, as shown by
$\widetilde{\mathbf{\bm{a}}}-\mathbf{\bm{b}}_{a}=\mathbf{\bm{n}}_{a}$ from Eq.
(1) and
$\widetilde{\mathbf{\bm{\omega}}}-\mathbf{\bm{b}}_{g}=\mathbf{\bm{n}}_{g}$
from Eq. (2), leading to accelerometer and gyroscope covariance matrices
$\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{a}}}}=({\sigma_{a}}^{2}/\Delta
t)\cdot\mathbf{\bm{I}}_{3\times 3}$ and
$\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{\omega}}}}=({\sigma_{g}}^{2}/\Delta
t)\cdot\mathbf{\bm{I}}_{3\times 3}$, respectively. Additional uncertainties
from ${}^{I_{0}}\widetilde{\mathbf{\bm{\omega}}}_{k}^{2}$ and
${}^{I_{0}}{\mathbf{\bm{\alpha}}}_{k}$ introduced by
${}^{I_{0}}\widehat{\mathbf{\bm{a}}}_{k}$ lead to covariance matrices
$\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{\omega}}}^{2}}=2\cdot\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{\omega}}}}^{2}$
and
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\alpha}}}=\sigma_{\mathbf{\bm{\alpha}}}^{2}\cdot\mathbf{\bm{I}}_{3\times
3}$, respectively. Consequently, the covariance matrix for
$\mathbf{\bm{r}}_{a}$ is
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{a}}=2\cdot\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{a}}}}+2\cdot\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{\omega}}}}^{2}+\mathbf{\bm{\Sigma}}_{\mathbf{\bm{\alpha}}}$
and that for $\mathbf{\bm{r}}_{g}$ is
$\mathbf{\bm{r}}_{g}=2\cdot\mathbf{\bm{\Sigma}}_{\widetilde{\mathbf{\bm{\omega}}}}$.
Residuals $\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}$ and
$\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}$ account for bias evolution:

$\displaystyle\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}$
$\displaystyle=\mathbf{\bm{b}}_{a_{n},k+1}-\mathbf{\bm{b}}_{a_{n},k}$
$\displaystyle\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}$
$\displaystyle=\mathbf{\bm{b}}_{g_{n},k+1}-\mathbf{\bm{b}}_{g_{n},k}.$

The covariance matrices for these residuals are
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{b}}_{a}}={\sigma_{\mathbf{\bm{b}}_{a}}}^{2}\Delta
t\cdot\mathbf{\bm{I}}_{3\times 3}$ and
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{b}}_{g}}={\sigma_{\mathbf{\bm{b}}_{g}}}^{2}\Delta
t\cdot\mathbf{\bm{I}}_{3\times 3}$, respectively.

###  IV-F Solution Approach

Upon identifying high-utility measurements, we formulate a nonlinear least-
square problem using the chosen measurement subset. Beyond the previously
mentioned residuals in Section IV-E, we incorporate additional residuals,
$\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}^{\ast}$ and
$\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}^{\ast}$, to bridge temporal gaps
between adjacent selected segments, denoted as
$\mathcal{D}^{info}=\\{\mathcal{D}^{m}|m\in\mathcal{M}\\}$ with
$\mathcal{M}=\\{\mathcal{M}_{1},\dotsc,\mathcal{M}_{M}\\}\subseteq\\{1,\dotsc,L\\}$:

$\displaystyle\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}^{\ast}$
$\displaystyle=\mathbf{\bm{b}}_{a_{n},(\mathcal{M}_{m+1}-1)\cdot
K+1}-\mathbf{\bm{b}}_{a_{n},\mathcal{M}_{m}\cdot K}$
$\displaystyle\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}^{\ast}$
$\displaystyle=\mathbf{\bm{b}}_{g_{n},(\mathcal{M}_{m+1}-1)\cdot
K+1}-\mathbf{\bm{b}}_{g_{n},\mathcal{M}_{m}\cdot K}.$

These additional residuals aim to maintain bias consistency across segments,
with their covariance matrices
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{b}}_{a}}^{\ast}$ and
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{b}}_{g}}^{\ast}$ modeled after a random
walk in discrete time as earlier discussed, resulting in
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{b}}_{a}}^{\ast}={\sigma_{\mathbf{\bm{b}}_{a}}}^{2}\cdot(\mathcal{M}_{m+1}-\mathcal{M}_{m})K\cdot\Delta
t\cdot\mathbf{\bm{I}}_{3\times 3}$ and
$\mathbf{\bm{\Sigma}}_{\mathbf{\bm{b}}_{g}}^{\ast}={\sigma_{\mathbf{\bm{b}}_{g}}}^{2}\cdot(\mathcal{M}_{m+1}-\mathcal{M}_{m})K\cdot\Delta
t\cdot\mathbf{\bm{I}}_{3\times 3}$.

These residuals allow for parameter estimation through solving a nonlinear
least-squares problem

$\begin{split}\min\Bigg{\\{}\sum_{l\in\mathcal{M}}&\Big{(}\sum_{\begin{subarray}{c}n\in\\{1,\dotsc,N\\}\\\
k\in\\{1,\dotsc,K\\}\end{subarray}}\big{(}\left\|\mathbf{\bm{r}}_{a}\right\|^{2}_{\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{a}}}+\left\|\mathbf{\bm{r}}_{g}\right\|^{2}_{\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{g}}}\big{)}\\\
+&\sum_{\begin{subarray}{c}n\in\\{0,\dotsc,N\\}\\\
k\in\\{1,\dotsc,K-1\\}\end{subarray}}\big{(}\left\|\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}\right\|^{2}_{\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}}}+\left\|\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}\right\|^{2}_{\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}}}\big{)}\Big{)}\\\
+\sum_{\begin{subarray}{c}m\in\\{1,\dotsc,M-1\\}\\\
n\in\\{0,\dotsc,N\\}\end{subarray}}&\big{(}\left\|\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}^{\ast}\right\|^{2}_{\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{a}}}^{\ast}}+\left\|\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}^{\ast}\right\|^{2}_{\mathbf{\bm{\Sigma}}_{\mathbf{\bm{r}}_{\mathbf{\bm{b}}_{g}}}^{\ast}}\big{)}\Bigg{\\}},\end{split}$

(3)

where $\left\|\cdot\right\|^{2}_{\mathbf{\bm{\Sigma}}}$ represents the
Mahalanobis distance using covariance matrix $\mathbf{\bm{\Sigma}}$.

##  V Simulation Experiments

In this section, we assess multi-IMU extrinsic self-calibration using
simulated trajectories from OpenVINS [20]. First, in Section V-B, we present
an edge case highlighting the need for adopting the greedy algorithm for
subset selection over M-largest selection, also noting the original greedy
algorithm (Alg. 1)’s efficiency issue. A sensitivity analysis in Section V-C
examines the impact of changes in parameter estimates on the measure of
information—a function of the parameter estimates. This analysis supports the
use of Alg. 2 as an alternative to Alg. 1. Section V-D compares the
calibration performances of Alg. 1 and Alg. 2 against the baseline of
calibrating with the full set of measurements, showing Alg. 2 significantly
reduces runtime—from minutes to roughly a quarter minute—without sacrificing
accuracy.

###  V-A Implementation Details

TABLE III: Reference IMU poses in simulation

index | position [cm] | orientation [deg]*
---|---|---
IMU0 | [0, 0, 0] | [0, 0, 0]
IMU1 | [20, 0, 0] | [180, 0, 0]
IMU2 | [0, 20, 0] | [0, 180, 0]
IMU3 | [0, 0, 20] | [0, 0, 180]
|  |

TABLE IV: *

* in this table, orientation is given as XYZ Euler angles

Table IV shows the reference poses for four IMUs mounted on a rigid body.
Accelerometer noise and bias instability were set at $\sigma_{a}=2\times
10^{-3}~{}\textrm{m}/\textrm{s}^{2}/\sqrt{\textrm{Hz}}$ and
$\sigma_{\mathbf{\bm{b}}_{a}}=3\times
10^{-3}~{}\textrm{m}/\textrm{s}^{2}\cdot\sqrt{\textrm{Hz}}$, respectively,
while for gyroscopes, these values were $\sigma_{g}=1.6968\times
10^{-4}~{}\textrm{rad}/\textrm{s}/\sqrt{\textrm{Hz}}$ and
$\sigma_{\mathbf{\bm{b}}_{g}}=1.9393\times
10^{-5}~{}\textrm{rad}/\textrm{s}\cdot\sqrt{\textrm{Hz}}$. Gyroscope
misalignment was generated by rotating each IMU’s reference orientation about
a uniformly random axis by an angle from a zero-mean normal distribution with
a standard deviation of $1^{\circ}$. IMU measurements were generated at 100
Hz.

###  V-B Limitations of Existing Methods for Self-Calibration with Informative
Subset Selection

TABLE V: Error in estimated extrinsic parameters and runtime for baseline,
M-largest selection, and greedy algorithm

Baseline | M-largest Selection |  Greedy Algorithm
---|---|---
$\mathbf{\bm{p}}$ [cm] |  $\mathbf{\bm{q}}$ [deg] |  $\mathbf{\bm{q}}_{g}$ [deg] |  $t$ [s] |  $\mathbf{\bm{p}}$ [cm] |  $\mathbf{\bm{q}}$ [deg] |  $\mathbf{\bm{q}}_{g}$ [deg] |  $t$ [s] |  $\mathbf{\bm{p}}$ [cm] |  $\mathbf{\bm{q}}$ [deg] |  $\mathbf{\bm{q}}_{g}$ [deg] |  $t$ [s]
0.0151 | 0.3685 | 2.1002 | 7.58 | 19.8605 | 1.4150 | 62.3996 | 1.91 | 0.0229 | 0.3090 | 0.8152 | 8.04
|  |  |  |  |  |  |  |  |  |  |

Figure 1: The intervals of subset within a roll-pitch-yaw profile, which is
used by the calibration process for the baseline (left), M-largest selection
(middle), and greedy algorithm (right), are highlighted.

We assess multi-IMU extrinsic self-calibration using three measurement
subsets: the full set (“baseline”), and subsets selected by M-largest
selection and the greedy algorithm. We simulated a trajectory with sinusoidal
oscillations in XYZ Euler angles every 20 seconds, maintaining consistent
information for each interval:

  * •

For $0\leq t<20$ seconds: $[\text{sin}(1\pi t),0,0]~{}\textrm{[rad]}$.

  * •

For $20\leq t<40$ seconds: $[0,\text{sin}(2\pi t),0]~{}\textrm{[rad]}$.

  * •

For $40\leq t<60$ seconds: $[0,0,\text{sin}(3\pi t)]~{}\textrm{[rad]}$.

For the greedy algorithm, the utility threshold $\lambda$ was set to 0.5,
following previous studies [14, 15]. For M-largest selection, we chose 15
segments, ensuring a similar subset size to the greedy algorithm. Initial
parameter estimates $\mathbf{\bm{p}}_{{I}_{n}}$, $\mathbf{\bm{q}}_{{I}_{n}}$,
and $\mathbf{\bm{q}}_{g_{n}}$ were set to zero, with angular acceleration
derived from numerical differentiation of IMU 0’s gyroscope measurements and
time-varying biases starting at zero.

Table V presents the calibration error and runtime for each process, and Fig.
1 illustrates the chosen segments. M-largest selection shows significantly
higher calibration errors than both the baseline and greedy algorithm. While
greedy algorithm selects new segments with utility (notably around 20s and
40s), the failure of M-largest selection in this regard likely results in
poorer estimates. Despite lower errors, greedy algorithm’s runtime is longer
than M-largest selection and the baseline, highlighting the need for the
efficient subset selection as introduced in Section III.

###  V-C Sensitivity Analysis of Fisher Information Matrix

The underlying assumption of Alg. 2 is the Fisher information matrix
($\mathbf{\bm{\mathcal{I}}}=\mathbf{\bm{J}}^{T}|_{\mathbf{\bm{\theta}}}\mathbf{\bm{J}}|_{\mathbf{\bm{\theta}}}$)
— a function of the estimate for the parameters $\mathbf{\bm{\theta}}$ — is
insensitive to the parameter estimate. To investigate this, we examined the
Spearman’s correlation coefficient [21] between the distribution of
information
$-1/2\cdot\log\big{|}\big{[}{\mathbf{\bm{\mathcal{I}}}}^{-1}\big{]}_{\mathbf{\bm{\Theta}}\mathbf{\bm{\Theta}}}\big{|}$
at reference extrinsic parameters (original distribution) and at values
altered by deviations in either position ($\delta{\mathbf{\bm{p}}}$) or
orientation ($\delta{\mathbf{\bm{q}}}$) (perturbed distribution). The study
utilized simulations across 48 different trajectories from five datasets
provided by OpenVINS [20].

Fig. 2 shows changes in the correlation coefficient are smaller for positional
deviations ($\delta{\mathbf{\bm{p}}}$) but larger for orientational deviations
($\delta{\mathbf{\bm{q}}}$), even with a few degrees. This indicates
information sensitivity mainly lies with orientation guesses rather than
position. Hence, ensuring accurate initial orientation guesses, as detailed in
Section IV-D, is key for effective segment selection by Alg. 2.

Figure 2: The correlation coefficient of the information as a function of the
difference ($\delta\mathbf{\bm{p}}$, $\delta\mathbf{\bm{q}}$) in guessed and
reference extrinsic parameters in the multi-IMU system. A coefficient of 1
indicates perfect correlation with the true value, while 0 denotes
independence.

###  V-D Comparison of Self-Calibration with Greedy Algorithm against
Benchmarks

TABLE VI: Ratio of selected segments and runtime for baseline, original
greedy, and greedy with utility evaluation at initial calibration parameters
over each dataset

| Ratio of selected segments [$\%$] |  Runtime [s]
---|---|---
Dataset |  $\\#$ of traj. | Length [m] | Duration [s] | Baseline | Greedy (Original) | Greedy (Init-Param) | Baseline | Greedy (Original) | Greedy (Init-Param)
uzh_fpv | 15 | 210.02 $\pm$ 144.69 | 39.03 $\pm$ 17.68 | 100.00 $\pm$ 0.00 | 59.30 $\pm$ 17.37 | 66.15 $\pm$ 21.01 | 5.98 $\pm$ 3.46 | 7.12 $\pm$ 4.44 | 3.41 $\pm$ 0.94
euroc_mav | 11 | 81.22 $\pm$ 23.44 | 122.58 $\pm$ 27.58 | 100.00 $\pm$ 0.00 | 26.06 $\pm$ 8.43 | 24.88 $\pm$ 7.41 | 18.21 $\pm$ 4.78 | 23.16 $\pm$ 6.22 | 5.72 $\pm$ 1.16
tum_vi | 6 | 115.25 $\pm$ 36.98 | 135.10 $\pm$ 12.52 | 100.00 $\pm$ 0.00 | 25.94 $\pm$ 7.03 | 26.33 $\pm$ 5.11 | 19.82 $\pm$ 1.99 | 32.45 $\pm$ 10.06 | 5.72 $\pm$ 1.06
kaist_vio | 11 | 37.10 $\pm$ 16.54 | 179.44 $\pm$ 30.00 | 100.00 $\pm$ 0.00 | 17.76 $\pm$ 9.74 | 17.66 $\pm$ 6.39 | 25.28 $\pm$ 11.88 | 30.91 $\pm$ 19.29 | 5.32 $\pm$ 1.79
kaist | 5 | 9618.87 $\pm$ 1958.50 | 1533.96 $\pm$ 656.45 | 100.00 $\pm$ 0.00 | 3.87 $\pm$ 2.26 | 3.44 $\pm$ 2.11 | 379.42 $\pm$ 135.80 | 460.81 $\pm$ 211.07 | 20.08 $\pm$ 6.60
|  |  |  |  |  |  |  |  |

TABLE VII: Runtime breakdown for baseline, original greedy, and greedy with
utility evaluation at initial calibration parameters over each dataset

| Baseline | Greedy (Original) | Greedy (Init-Param)
---|---|---|---
Dataset | Evaluate | Calibrate | Total | Evaluate | Calibrate | Total | Evaluate | Calibrate | Total
uzh_fpv | - | 5.94 $\pm$ 3.43 | 5.98 $\pm$ 3.46 | 2.52 $\pm$ 2.05 | 4.50 $\pm$ 2.39 | 7.12 $\pm$ 4.44 | 0.15 $\pm$ 0.07 | 3.19 $\pm$ 0.86 | 3.41 $\pm$ 0.94
euroc_mav | - | 18.09 $\pm$ 4.76 | 18.21 $\pm$ 4.78 | 13.84 $\pm$ 4.11 | 8.91 $\pm$ 2.94 | 23.16 $\pm$ 6.22 | 0.46 $\pm$ 0.13 | 5.12 $\pm$ 1.19 | 5.72 $\pm$ 1.16
tum_vi | - | 19.69 $\pm$ 1.98 | 19.82 $\pm$ 1.99 | 18.50 $\pm$ 6.46 | 13.48 $\pm$ 3.89 | 32.45 $\pm$ 10.06 | 0.54 $\pm$ 0.06 | 5.04 $\pm$ 1.06 | 5.72 $\pm$ 1.06
kaist_vio | - | 25.13 $\pm$ 11.82 | 25.28 $\pm$ 11.88 | 21.11 $\pm$ 14.18 | 9.27 $\pm$ 5.43 | 30.91 $\pm$ 19.29 | 0.58 $\pm$ 0.26 | 4.56 $\pm$ 1.51 | 5.32 $\pm$ 1.79
kaist | - | 377.85 $\pm$ 135.22 | 379.42 $\pm$ 135.80 | 371.62 $\pm$ 172.61 | 81.47 $\pm$ 41.68 | 460.81 $\pm$ 211.07 | 4.61 $\pm$ 1.98 | 13.46 $\pm$ 5.33 | 20.08 $\pm$ 6.60
|  |  |  |  |  |  |  |  |

TABLE VIII: Memory footprint (RAM usage in MB) for baseline, original greedy,
and greedy with utility evaluation at initial calibration parameters over each
dataset

Dataset | Baseline | Greedy (Original) | Greedy (Init-Param)
---|---|---|---
uzh_fpv | 92.69 $\pm$ 40.62 | 55.38 $\pm$ 14.88 | 65.13 $\pm$ 32.34
euroc_mav | 235.46 $\pm$ 50.42 | 82.70 $\pm$ 14.43 | 69.10 $\pm$ 33.19
tum_vi | 264.21 $\pm$ 24.86 | 92.54 $\pm$ 17.72 | 100.73 $\pm$ 38.73
kaist_vio | 274.96 $\pm$ 124.70 | 80.69 $\pm$ 45.38 | 85.89 $\pm$ 47.05
kaist | 1678.32 $\pm$ 350.95 | 141.08 $\pm$ 11.22 | 135.75 $\pm$ 8.31
|  |  |

Figure 3: Absolute error in estimated extrinsic parameters for baseline,
original greedy, and greedy with utility evaluation at initial calibration
parameters over each dataset in simulation.

We assess the self-calibration with subset selection in Alg. 2 (“Greedy (Init-
Param)”) against Alg. 1 (“Greedy (Original)”) and a “Baseline” using all
measurements. Our focus lies on runtime and calibration performance
comparisons across 48 trajectories within five OpenVINS-supported datasets.
For both the greedy methods, information evaluation occurred every second,
using a $\lambda=0.5$ threshold, with initial parameter guesses aligned with
the earlier discussion.

Table VI presents the chosen measurement subset ratio and runtime for each
calibration process, with Fig. 3 showing the absolute extrinsic parameter
errors for all IMUs. The runtime breakdown and memory usage are presented in
Tables VII and VIII, respectively. Greedy (Init-Param) consistently reduced
runtime, with significant improvements in the kaist dataset, while Greedy
(Original) showed less or negative improvement. Calibration accuracy remained
consistent both for greedy methods and the baseline, achieving sub-centimeter
and sub-degree precision, except in the kaist dataset. Greedy methods
particularly outperformed the baseline in the kaist, which was collected in a
driving car and primarily consists of low-utility measurements, by excluding
these low-utility measurements to improve problem observability. The subset
selection addresses the challenges of unobservable or weakly observable
calibration scenarios that can arise in self-calibration with full
measurements.

##  VI Hardware Experiments

This section evaluates self-calibration with measurement subsets using real-
world data from a dataset featured in our prior work [13]. As in simulations,
a sensitivity analysis in Section VI-B examines how changes in parameter
estimates affect the measure of information, supporting the use of Alg. 2 as
an alternative to Alg. 1. Section VI-C compares the calibration performances
of Alg. 1 and Alg. 2 with a full-measurement baseline, highlighting Alg. 2’s
reduction in runtime—from approximately one minute to just a second—without
loosing calibration accuracy.

###  VI-A Implementation Details

Figure 4: Sensor rig used for data collection.

Fig. 4 shows the sensor rig used for data collection, featuring two IMUs (IMU
0 and IMU 1) and a camera (CAM), with only the IMU data utilized for
evaluation. The reference relative position $\mathbf{\bm{p}}_{I_{1}}$ and
orientation $\mathbf{\bm{q}}_{I_{1}}$ between the IMUs are
$[-10,-10,0]\;[\textrm{cm}]$ and $[0,0,0]\;[\textrm{deg}]$ in XYZ Euler
angles, with assumed zero gyroscope misalignment; $\mathbf{\bm{q}}_{g_{0}}$
and $\mathbf{\bm{q}}_{g_{1}}$ are also $[0,0,0]\;[\textrm{deg}]$. Allan
variance analysis provided white noise and bias instability of the
accelerometers $\sigma_{a}=3.3\times
10^{-3}~{}\textrm{m}/\textrm{s}^{2}/\sqrt{\textrm{Hz}}$ and
$\sigma_{\mathbf{\bm{b}}_{a}}=2.4\times
10^{-4}~{}\textrm{m}/\textrm{s}^{2}\cdot\sqrt{\textrm{Hz}}$ and of the
gyroscopes $\sigma_{g}=8.3\times
10^{-5}~{}\textrm{rad}/\textrm{s}/\sqrt{\textrm{Hz}}$ and
$\sigma_{\mathbf{\bm{b}}_{g}}=3.7\times
10^{-6}~{}\textrm{rad}/\textrm{s}\cdot\sqrt{\textrm{Hz}}$.

Acknowledging that comparing estimated calibration parameters to reference
values might not entirely reflect accuracy, as noted in our previous work, we
also measure reprojection error between the IMUs, with a lower error
indicating more precise estimation.

###  VI-B Sensitivity Analysis of Fisher Information Matrix

As in the sensitivity analysis conducted in simulations (in Section V-C), we
assessed how deviations from true calibration parameters affect the measured
information using the hardware data. In Fig.5, we examined the Spearman’s
correlation coefficient [21] between the information distribution,
$-1/2\cdot\log\big{|}\big{[}{\mathbf{\bm{\mathcal{I}}}}^{-1}\big{]}_{\mathbf{\bm{\Theta}}\mathbf{\bm{\Theta}}}\big{|}$,
evaluated at reference extrinsic parameters and at deviations in either
relative position ($\delta{\mathbf{\bm{p}}}$) or orientation
($\delta{\mathbf{\bm{q}}}$).

Similar to that in simulations, sensitivity to orientation deviations
($\delta{\mathbf{\bm{q}}}$) was evident, showing a greater decline than for
position deviations ($\delta{\mathbf{\bm{p}}}$). This highlights the
importance of good initial guesses for relative orientation, as advocated in
Section IV-D, for sucessful measurement subset selection.

Figure 5: The correlation coefficient of the information as a function of the
difference ($\delta\mathbf{\bm{p}}$, $\delta\mathbf{\bm{q}}$) in guessed and
reference extrinsic parameters in the two-IMU system on hardware. A
coefficient of 1 indicates perfect correlation with the true value, while 0
denotes independence.

###  VI-C Comparison of Self-Calibration with Efficient Greedy Algorithm
against Benchmarks

TABLE IX: Estimated extrinsics, segment selection ratio, and runtime for
baseline, original greedy, and greedy with utility evaluation at initial
calibration parameters on hardware data

Trajectory |  | Baseline | Greedy (Original) | Greedy (Init-Param)
---|---|---|---|---
baseline (1274 [s]) |  $\mathbf{\bm{p}}_{I_{1}}$ [cm] | [-9.80 -9.82 -0.24] | [-9.85 -9.63 -0.34] | [-9.64 -9.58 -0.29]
$\mathbf{\bm{q}}_{I_{1}}$ [deg] | [4.15 3.28 0.27] | [ 3.39 2.01 -0.42] | [ 1.98 0.61 -0.40]
$\mathbf{\bm{q}}_{g_{0}}$ [deg] | [-0.74 -0.45 0.04] | [ 2.47 4.34 -1.08] | [-0.90 1.06 0.21]
$\mathbf{\bm{q}}_{g_{1}}$ [deg] | [3.11 3.31 0.45] | [ 5.62 6.82 -1.05] | [0.84 2.12 0.09]
Selected segments [%] | 100.00 | 0.73 | 1.14
Runtime [s] | 64.09 | 16.15 | 0.98
blurry (1388 [s]) |  $\mathbf{\bm{p}}_{I_{1}}$ [cm] | [-9.81 -9.77 -0.26] | [-9.84 -9.78 -0.30] | [-9.79 -9.72 -0.27]
$\mathbf{\bm{q}}_{I_{1}}$ [deg] | [5.27 4.27 0.25] | [ 2.22 0.98 -0.25] | [ 1.96 0.96 -0.65]
$\mathbf{\bm{q}}_{g_{0}}$ [deg] | [-2.96 -2.74 0.00] | [-0.43 -0.06 -0.27] | [-1.76 -1.13 0.45]
$\mathbf{\bm{q}}_{g_{1}}$ [deg] | [2.03 2.11 0.20] | [ 1.65 1.28 -0.19] | [0.09 0.15 0.05]
Selected segments [%] | 100.00 | 1.94 | 0.82
Runtime [s] | 75.81 | 51.01 | 0.90
ill-lit (1276 [s]) |  $\mathbf{\bm{p}}_{I_{1}}$ [cm] | [-9.82 -9.86 -0.21] | [-9.67 -9.96 -0.14] | [-9.65 -9.83 -0.17]
$\mathbf{\bm{q}}_{I_{1}}$ [deg] | [-5.03 -5.79 0.28] | [ 2.50 1.29 -0.02] | [5.03 3.50 0.10]
$\mathbf{\bm{q}}_{g_{0}}$ [deg] | [-0.36 -0.34 0.02] | [ 1.70 0.69 -0.18] | [ 2.55 1.89 -0.40]
$\mathbf{\bm{q}}_{g_{1}}$ [deg] | [-5.79 -5.58 0.39] | [3.82 2.42 0.02] | [7.24 6.20 0.21]
Selected segments [%] | 100.00 | 2.19 | 1.05
Runtime [s] | 73.75 | 59.37 | 0.99
|  |  |  |

TABLE X: IMU measurement reprojection errors for baseline, original greedy,
and greedy with utility evaluation at initial calibration parameters on
hardware data

| Dataset | Baseline | Greedy (Original) | Greedy (Init-Param)
---|---|---|---|---
accelerometer ([$\textrm{m}/\textrm{s}^{2}$]) | baseline | 0.8160 $\pm$ 0.4947 | 0.6015 $\pm$ 0.4947 | 0.4749 $\pm$ 0.4760
blurry | 1.0981 $\pm$ 0.5470 | 0.5069 $\pm$ 0.5470 | 0.5041 $\pm$ 0.4257
ill-lit | 1.6424 $\pm$ 0.3951 | 0.4697 $\pm$ 0.3951 | 0.8968 $\pm$ 0.3911
gyroscope ([$\textrm{rad}/\textrm{s}$]) | baseline | 0.0942 $\pm$ 0.0139 | 0.0942 $\pm$ 0.0139 | 0.0942 $\pm$ 0.0141
blurry | 0.0996 $\pm$ 0.0209 | 0.0996 $\pm$ 0.0219 | 0.0996 $\pm$ 0.0219
ill-lit | 0.0926 $\pm$ 0.0101 | 0.0925 $\pm$ 0.0103 | 0.0930 $\pm$ 0.0107
|  |  |  |

Table IX shows estimated calibration parameters $\mathbf{\bm{p}}_{I_{1}}$,
$\mathbf{\bm{q}}_{I_{1}}$, $\mathbf{\bm{q}}_{g_{0}}$, and
$\mathbf{\bm{q}}_{g_{1}}$ along with the chosen measurement subset ratio and
runtime for each calibration process. Table X shows the IMU reprojection error
from the three methods across trajectories, with significant differences (via
t-test) in bold for the greedy methods. Initial parameter guesses were as
previously discussed in simulations. The results indicate the greedy
algorithms use less than 3% of total measurements to achieve calibration
results that closely align with reference values. Greedy (Init-Param)
significantly reduces runtime from over a minute to approximately a second,
whereas Greedy (Original) achieves a modest reduction, cutting the runtime by
about fifteen seconds from the baseline.

##  VII Conclusion

This paper presents a multi-IMU extrinsic calibration by selecting high-
utility measurement subsets. We show that, within our systems of interest,
utility—a function of the parameter estimates—is insensitive to the changes in
their values, thus allowing its assessment at some initial guess. This
eliminates the need for frequent recalibrations during subset selection,
significantly reducing calibration computation time. Our simulations and
experiments confirm this efficiency gain, suggesting the method’s
applicability for self-calibration on resource-limited platforms where
handling large datasets efficiently is crucial.

##  Appendix A Proof of Time Complexity of Original Greedy Algorithm

In this appendix, we derive the time complexity of the original greedy
algorithm (Alg. 1).

We will use the following symbols:

  * •

$L$ means the number of total segments.

  * •

$M$ means the number of finally chosen segments (i.e., $M\leq L$).

  * •

$K$ means the number of timesteps (i.e., data) in each segment.

  * •

$p$ means the number of parameters being estimated.

We assume $K>p$ in general.

  1. 1.

Line 5 in Alg. 1 iterates $L$ times.

  2. 2.

Line 6 conducts nonlinear least-square on $m\times K$ data, where $m$ is the
number of segments in $\mathcal{D}^{info}\cup\mathcal{D}^{new}$
($\mathcal{O}(mK\cdot p^{2})$).

  3. 3.

Line 7 evaluates Jacobians
$\mathbf{\bm{J}}(\mathcal{D}^{info})_{\mathbf{\bm{\hat{\theta}}}^{+}}$ (size:
$(m-1)K\times p$) and
$\mathbf{\bm{J}}(\mathcal{D}^{new})_{\mathbf{\bm{\hat{\theta}}}^{+}}$ (size:
$K\times p$) and computes the Fisher information matrix
($\mathcal{O}(mKp^{2})$).

  4. 4.

Lines 8-9 calculate the Fisher information matrix’s inverse and determinant
(both $\mathcal{O}(p^{3})$).

Operations on lines 6 and 7 within the loop are computational bottlenecks. The
sum of $m$ in the loop can reach $L^{2}$ in the worst case and $L$ in the best
case, leading to Alg. 1’s time complexity expressed as:

$T(n)=\mathcal{O}(\sum_{m}mKp^{2})=\begin{cases}\mathcal{O}(LKp^{2})&\textrm{in
the best}\\\ \mathcal{O}(L^{2}Kp^{2})&\textrm{in the worst}\\\ \end{cases}$

##  Appendix B Proof of Time Complexity of Greedy Algorithm Using Only Initial
Calibration Parameters

In this appendix, we derive the time complexity of the greedy algorithm with
utility evaluation at initial calibration parameters (Alg. 2). We use the same
symbols as introduced in Appendix A.

  1. 1.

Line 4 in Alg. 2 iterates $L$ times.

  2. 2.

Line 5 evaluates Jacobian
$\mathbf{\bm{J}}(\mathcal{D}^{new})|_{\mathbf{\bm{\theta}}^{0}}$ (size:
$K\times p$) and computes the Fisher information matrix
($\mathcal{O}(Kp^{2})$).

  3. 3.

Lines 6-7 calculates the Fisher information matrix’s inverse and determinant
(both $\mathcal{O}(p^{3})$).

  4. 4.

After the loop, line 12 conducts nonlinear least-square on $M\times K$ data
($\mathcal{O}(MK\cdot p^{2})$).

The operation on line 5 within the loop is a computational bottleneck, leading
to Alg. 2’s time complexity expressed as:

$T(n)=\mathcal{O}(L\times Kp^{2}+(MK)p^{2})=\mathcal{O}(LKp^{2}).$

## References

  * [1] J.-O. Nilsson and I. Skog, “Inertial sensor arrays—a literature review,” in _Proc. IEEE Eur. Navigation Conf_ , 2016.
  * [2] I. Skog, J.-O. Nilsson, P. Händel, and A. Nehorai, “Inertial sensor arrays, maximum likelihood, and cramér–rao bound,” _IEEE Trans. Signal Process._ , vol. 64, no. 16, pp. 4218–4227, 2016.
  * [3] K. Parsa, J. Angeles, and A. K. Misra, “Estimation of the flexural states of a macro-micro manipulator using point-acceleration data,” _IEEE Trans. Robot._ , vol. 21, no. 4, pp. 565–573, 2005.
  * [4] I. B. Wijayasinghe, M. N. Saadatzi, S. Abubakar, and D. O. Popa, “A study on optimal placement of accelerometers for pose estimation of a robot arm,” in _Proc. IEEE Int. Conf. Robot. Automat._ , 2018, pp. 1444–1451.
  * [5] M. Zhang, X. Xu, Y. Chen, and M. Li, “A lightweight and accurate localization algorithm using multiple inertial measurement units,” _IEEE Robot. Autom. Lett._ , vol. 5, no. 2, pp. 1508–1515, 2020.
  * [6] K. Eckenhoff, P. Geneva, and G. Huang, “Mimc-vins: A versatile and resilient multi-imu multi-camera visual-inertial navigation system,” _IEEE Trans. Robot._ , vol. 37, no. 5, pp. 1360–1380, 2021.
  * [7] S. Y. Cho and C. G. Park, “A calibration technique for a redundant imu containing low-grade inertial sensors,” _ETRI J._ , vol. 27, no. 4, pp. 418–426, 2005.
  * [8] P. Schopp, L. Klingbeil, C. Peters, and Y. Manoli, “Design, geometry evaluation, and calibration of a gyroscope-free inertial measurement unit,” _Sensors Actuators A: Phys._ , vol. 162, no. 2, pp. 379–387, 2010.
  * [9] K. He, J. Han, and Y. Shao, “A novel redundant inertial measurement unit and calibration algorithm,” in _Proc. IEEE Int. Conf. Optoelectron. Microelectronics_ , 2013, pp. 18–23.
  * [10] J. Rehder, J. Nikolic, T. Schneider, T. Hinzmann, and R. Siegwart, “Extending kalibr: Calibrating the extrinsics of multiple imus and of individual axes,” in _Proc. IEEE Int. Conf. Robot. Automat._ , 2016, pp. 4304–4311.
  * [11] P. Schopp, H. Graf, W. Burgard, and Y. Manoli, “Self-calibration of accelerometer arrays,” _IEEE Trans. Instrum. Meas._ , vol. 65, no. 8, pp. 1913–1925, 2016.
  * [12] D. Kim, S. Shin, and I. S. Kweon, “On-line initialization and extrinsic calibration of an inertial navigation system with a relative preintegration method on manifold,” _IEEE Trans. Autom. Sci. Eng._ , vol. 15, no. 3, pp. 1272–1285, 2017.
  * [13] J. Lee, D. Hanley, and T. Bretl, “Extrinsic calibration of multiple inertial sensors from arbitrary trajectories,” _IEEE Robot. Autom. Lett._ , vol. 7, no. 2, pp. 2055–2062, 2022.
  * [14] J. Maye, P. Furgale, and R. Siegwart, “Self-supervised calibration for robotic systems,” in _Proc. IEEE Intell. Veh. Symp._ , 2013, pp. 473–480.
  * [15] J. Maye, H. Sommer, G. Agamennoni, R. Siegwart, and P. Furgale, “Online self-calibration for robotic systems,” _Int. J. Robot. Res_ , vol. 35, no. 4, pp. 357–380, 2016.
  * [16] C. R. Rao, “Minimum variance and the estimation of several parameters,” in _Math. Proc. Camb. Phil. Soc._ , vol. 43, no. 2. Cambridge University Press, 1947, pp. 280–283.
  * [17] T. Schneider, M. Li, C. Cadena, J. Nieto, and R. Siegwart, “Observability-aware self-calibration of visual and inertial sensors for ego-motion estimation,” _IEEE Sensors J._ , vol. 19, no. 10, pp. 3846–3860, 2019.
  * [18] J. Lv, X. Zuo, K. Hu, J. Xu, G. Huang, and Y. Liu, “Observability-aware intrinsic and extrinsic calibration of lidar-imu systems,” _IEEE Trans. Robot._ , 2022.
  * [19] Z. Yang and S. Shen, “Monocular visual–inertial state estimation with online initialization and camera–imu extrinsic calibration,” _IEEE Trans. Autom. Sci. Eng._ , vol. 14, no. 1, pp. 39–51, 2016.
  * [20] P. Geneva, K. Eckenhoff, W. Lee, Y. Yang, and G. Huang, “Openvins: A research platform for visual-inertial estimation,” in _Proc. IEEE Int. Conf. Robot. Automat._ , 2020, pp. 4666–4672.
  * [21] J. H. Zar, “Spearman rank correlation: overview,” _Wiley StatsRef: Statistics Reference Online_ , 2014.