Introduction

Inertial Measurement Unit (IMU's) have a profound role in the control of mobile vehicles, robotic manipulators, capturing human body motions and augmented reality applications in smartphones. Fused with a magnetometer, an IMU enables precise orientation estimates. In addition, IMU's provide velocity and position change estimates over short time periods. The quality of those estimates largely depend on correct sensor models. Hence, estimating the model parameters, known as calibration, is an important procedure that has major effect on the eventual usability of the IMU. The upswing of micro-electromechanical (MEMS) based IMU's resulted in many scientific publications describing calibration methods [6,10]. Traditionally an external system is used that applies suitable reference signals for a one-time calibration procedure right after the IMU was manufactured.

However, environmental influences, like temperature changes and mechanical stress, cause deviations of the true parameter values over time. Hence, an easy to perform, without the need for extra equipment, calibration procedure is desired to correct for those recurred systematic errors.

Supported by the Dutch SIA RAAK Postdoc program

While most literature focusses on the calibration of a single sensor, or only the gyroscope and accelerometer triads, do some include the magnetometer as well [8]. However, the complementary relations between inertial and magnetic signals are often minimally exploited [2,4]. Estimating the sensors jointly, without the need to perform complex movements or equipment, was only recently published by Chow et.al. [3]. It illustrates the possibilities of a powerful, yet easy to use calibration procedure. However, their choice for marginalization the navigation states puts some restrictions on the prior, the movement, its bandwidth and outlier handling.

We propose a similar method, yet different approach for an easy in-use calibration of consumer grade triaxial inertial and magnetometer sensors without the need for external equipment. The novelty can be found in various aspects:

-A tight coupled, batched, sensor fusion approach is used to benefit from all sensor observations in a joint manner. -Flexible in the number of rigidly connected IMU's (MIMU), their sampling rates, and the motion trajectory. -Exploiting the centripletal and tangential forces of a MIMU constellation.

-Simultaneous estimation of all inertial and magnetic calibration parameters.

-Simultaneous estimation of the navigation states.

-Sample based handling of measurement outliers.

-Robust for any magnetic disturbance and temporally or spatially inhomogenous fields.

This paper is structured as follows: we will first outline the theoretical framework and highlight the method. Then, the method is demonstrated in three different situations. First, a MIMU array is calibrated, where the orientation state is compared with an optical reference system. In the second situation a MIMU array is calibrated using multiple static periods and compared with an existing popular calibration approach. Third, the method is applied on an IMU embodied in an Apple iPhone X during typical smartphone usage.

Method

This section elaborates on the sensor, motion and constraint models. Diagonal gain matrices are indicated with K, lower triangular matrices describing the misalignments between the senstive axes with N , biases vectors with b and Gaussian noise vectors with e [11]. Super or sub-script upper-case letters represent the local world frame L, module body frame B, magnetometer M i , gyroscope G i , sensor frame S i . Orientation are indicated with an orientation matrix R.

We assume a rigid body on which one or multiple sensors (S i ) are attached. Hence, the sensor's orientation with respect to the global frame can be written as:

R SiL t = R SiB R BL t (1)

Sensor models

Accelerometer: The output of an accelerometer at time t is modeled as the sum of the linear a Si Si,t and gravitational acceleration g Si :

y Si a,t = K Si a N Si a R SiB a B Si,t + R SiL t g L + b Si a + e a,t , e a,t ∼ N (0, Σ a ) (2)

This linear acceleration at pick-up point (S i ) is a summation of the body, centripetal and tangential accelerations. Latter can be expressed using the lever-arm p B Si , angular velocity and angular acceleration:

a B Si,t = R BL t a L B,t + α B LB,t × p B Si + ω B LB,t × ω B LB,t × p B Si (3)

where R GiSi is the relative orientation, or triad misalignment, between the accelerometer and gyroscope pair.

Gyroscope: The output of a gyroscope is modeled as an angular velocity measured in the gyroscope frame with respect to the inertial frame (ω Gi

LGi,t ). We assume that the gyroscope and accelerometer are rigidly connected to the underlying body, hence ω BSi = ω SiGi = 0:

y Gi g,t = K Gi g N Gi g R GiSi R SiB ω B LB,t + b Gi g + e g,t , e g,t ∼ N (0, Σ g )(4)

Magnetometer: The output of a magnetometer at time t is modeled as a local measured homogenous magnetic field m L which is affected by soft (D) and hard iron (o) effects:

y Mi m,t = D Mi R SiB R BL m L + o Mi m + e m,t , e m,t ∼ N (0, Σ m )(5)

It should be noted that the triad misalignment between magnetometer and accelerometer is captured in matrix D Mi and the magnetometer bias is captured in the offset vector o Mi . The local magnetic field is a function of the local magnetic dip angle θ and modeled as:

m L = cos θ 0 − sin θ T (6)

Motion and constraint models

Motion model: The kinematic relations of the position p, velocity v , acceleration a, orientation q, angular velocity ω and angular acceleration α are given by the following proces model which is driven by a zero mean acceleration noise model:

p L t+T = p L t + T v L t + T 2 2 a L t + w p v L t+T = v L t + T a L t + w v a L t+T = 0 + w a q LB t+T = q LB t exp T 2 ω B LB,t + T 2 α B LB,t + w q ω B LB,t+T = ω B LB,t + T α B LB,t + w ω α B LB,t+T = 0 + w α

where T is the sample period, is the quaternion product operator and exp the quaternion exponential. The process noises are being described by:

w X L B t ∼ N (0, Σ X ), X ∈ {p, v, a, q, ω, α, }(7)

It should be noted that w p , w v , w a represent the same linear acceleration uncertainty, and w q , w ω , w α the same angular acceleration uncertainty.

Zero net position:

The user is asked to return the IMU to approximately the same position as where it has been picked-up. This knowledge is modeled as a zero net position change measurement:

p t+Tp = p t + e np(8)

where T p is the time difference between pick-up and return.

Zero velocity: Whenever a zero movement period is detected using a movingvariance detector [9] on the inertial sensor readings, the angular and translational velocities are assumed to be zero:

v L t = 0 + e v0(9)ω B LB,t = 0 + e ω0 .(10)

Calibration procedure

The model parameters and navigation states are estimated simultaneously using a large scale iterative non-linear least squares solver [1]. By rewriting and summing the model terms, one can rephrase the entire cost function as a bundle adjustment problem. The following set of unknown are included:

1) Accelerometer and Gyroscope gains, misalignments, biases:

K Si a , K Gi g , N Si a , N Gi g , b Si a , b Gi g ∀i ∈ S(11)

2) Magnetometer gains, magnetometer offsets, and the magnetic dip angle:

D Mi , o Mi , θ ∀i ∈ S(12)

2) Boresights, lever arms and sensor triad misalignments:

R BSi , p B Si , R SiGi ∀i ∈ S(13)

2) Body kinematics:

p L B,t , v L B,t , a L B,t , q LB t , ω B LB,t , α B LB,t ∀t ∈ T (14)

where S denotes the set of IMU and magnetic sensors and T the set of time instances.

Outliers are efficiently detected and suppressed by embedding a Cauchy Loss function for each observation individually. Sensor covariances are either obtained from data-sheets or derived from an Allan Variance plot. Other covariance values were chosen such that they have a realistic meaning, e.g. sub-centimeter level zero net position change, millimeter/s and millidegree/s zero velocity uncertainties.

Experiments

Our approach has been evaluated in different indoor experiments. Before start of the recording sensors were running long enough to stabilize the temperature. In each experiment a different estimation aspect will be emphasized:

1. MIMU module, calibrated with two static poses. The orientation of the insample trial is compared with an optical reference. 2. MIMU module, calibrated with 20 static poses. The position of an out-ofsample trial is compared with a different calibration set and optical reference. 3. Smartphone IMU module, calibrated with two static periods, representing a typical smartphone movement. The pose of the in-sample trial is presented.

MIMU two static periods

An Inertial Elements MIMU22BT multi IMU, sampled at 62.5Hz, module was used, see Fig. 1. The MIMU was manipulated indoors in a rather arbitrary way for about one minute while its movements were recorded by a passive optical reference system (Vicon @100Hz). Processing the raw inertial and magnetic sensor readings yielded the estimated calibration parameters. Uncalibrated and calibrated sensor values are depicted in Fig. 2. In addition, the in-sample orientation estimate was compared with the optical reference, see Fig. 3.

MIMU 20 static periods

In our previous work we described a new method for Pedestrian Dead Reckoning (PDR) based on MIMU sensors [5]. For the experiments, a popular joint gyroacc calibration method described by Tedali et.al. [10] was used to calibrate to inertial sensors. We re-used both the calibration and experimental datasets. The calibration dataset was used to calibrate the sensor using the method described in this paper. The experimental data set contains the intertial data of a subject who traversed repeatedly an indoor oval path (800 m). Subsequently, the estimated trajectory using a Kalman based PDR approach [7] for the dataset that has been calibrated with the two differente parametersets is depicted in Fig. 4. Visible are an optical reference (left), an EKF estimate based on the original inertial calibration set (middle) and a reconstruction using a calibration set obtained from our method (right).

Smartphone, single IMU, 2 static periods

Typical smartphone usage includes a movement sequence of picking it up from a table, perform a phoning or texting activity, and finally place the phone back on the table top. We simulated this situation by picking up the phone, rotating it in different directions while translating in the global vertical-direction for about 30 cm, and finally return it to the initial resting spot. Figure 5a illustrates the in-sample estimates of the translational kinematics. In addition, in-sample orientation estimates are visible by the (un)calibrated magnetometer readings in Fig. 5b.

Conclusion and Discussion

This paper describes a novel method for the calibration of inertial and magnetic measurement units. It allows for simultaneous estimation of calibration parameter and navigation states and flexible in the number of sensors used.

The advantages, suitability and performance is demonstrated in different experiments. Especially experiment 2 shows the added value of proper estimated parameters when a MIMU module is used for PDR activities.

In a follow up study we would like to include the uncertainty on the parameter estimates by evaluating the Hessian matrix. In addition, the framework can be easily extended by inclusion of the sensor covariances and time-dependent bias states.