On Inertial Body Tracking in the Presence of Model Calibration Errors

Abstract

In inertial body tracking, the human body is commonly represented as a biomechanical model consisting of rigid segments with known lengths and connecting joints. The model state is then estimated via sensor fusion methods based on data from attached inertial measurement units (IMUs). This requires the relative poses of the IMUs w.r.t. the segments—the IMU-to-segment calibrations, subsequently called I2S calibrations—to be known. Since calibration methods based on static poses, movements and manual measurements are still the most widely used, potentially large human-induced calibration errors have to be expected. This work compares three newly developed/adapted extended Kalman filter (EKF) and optimization-based sensor fusion methods with an existing EKF-based method w.r.t. their segment orientation estimation accuracy in the presence of model calibration errors with and without using magnetometer information. While the existing EKF-based method uses a segment-centered kinematic chain biomechanical model and a constant angular acceleration motion model, the newly developed/adapted methods are all based on a free segments model, where each segment is represented with six degrees of freedom in the global frame. Moreover, these methods differ in the assumed motion model (constant angular acceleration, constant angular velocity, inertial data as control input), the state representation (segment-centered, IMU-centered) and the estimation method (EKF, sliding window optimization). In addition to the free segments representation, the optimization-based method also represents each IMU with six degrees of freedom in the global frame. In the evaluation on simulated and real data from a three segment model (an arm), the optimization-based method showed the smallest mean errors, standard deviations and maximum errors throughout all tests. It also showed the lowest dependency on magnetometer information and motion agility. Moreover, it was insensitive w.r.t. I2S position and segment length errors in the tested ranges. Errors in the I2S orientations were, however, linearly propagated into the estimated segment orientations. In the absence of magnetic disturbances, severe model calibration errors and fast motion changes, the newly developed IMU centered EKF-based method yielded comparable results with lower computational complexity.

1. Introduction

Inertial motion capturing has found widespread use in various applications, including biomechanics and health as two prominent ones [1,2,3,4]. This development is, among other reasons, driven by the availability of smaller, cheaper and more precise hardware [4]. Inertial measurement units (IMUs) comprise gyroscopes and accelerometers providing 3D acceleration and 3D rotational velocity. In most cases they also contain magnetometers adding 3D magnetic fields. This is also referred to as MIMUs. From these measurements motion information can be estimated through sensor fusion.

A considerable amount of literature deals with fusion techniques for orientation tracking based on a single IMU [5,6,7,8,9,10]. Especially the separation of body acceleration and acceleration due to gravity [6,8,9,11] and the handling of magnetic disturbances [6,10] received attention. Many applications, however, require knowledge about the motion of several (connected) body segments, often in terms of joint angles, as derived from multiple attached IMUs [12,13,14]. In this case, rather than independently tracking the orientation of each segment by means of the attached IMU, biomechanical models of the human body including IMU-to-segment placement are used in the sensor fusion method for reducing tracking errors [15,16,17,18,19,20]. Here, the body is typically modeled as a set of rigid segments (the bones) with known lengths, which are connected through frictionless joints of varying degrees of freedom (DoFs), and the IMUs are assumed to be rigidly mounted on the associated segments (mostly with a one-to-one mapping) [21]. The respective transformation, i.e., the orientation and position of each IMU w.r.t. the segment it is attached to, is called IMU-to-segment (subsequently abbreviated I2S) calibration [22].

It is obvious that such model assumptions introduce errors associated to these assumptions. In this work, we are focusing on the influence of selected model calibration errors, more specifically, on errors in the calibrated I2S orientations, positions and segment lengths, on the segment orientation estimation accuracy. Other calibration errors (e.g., joint rotation axes [23]) and errors associated to the simplification of the biomechanical model (e.g., limited joint DoFs [19] or soft tissue artifacts [24], which are motion reconstruction errors of the skeleton due to movement of e.g., skin or cloths w.r.t. the bone) are not in the focus here.

Sensor fusion algorithms were developed to better handle or compensate for sensor errors [8,9,20,25], such as noise and bias, and environmental effects, such as magnetic disturbances [6,10]. In contrast, model calibration errors were not intensively addressed in literature; though they represent a significant source of error [25] (cf. Section 1.2). Therefore, this work proposes three new/adapted sensor fusion methods with the aim of providing increased robustness against such model calibration errors. The performances of these methods, in combination with their dependence on magnetometer usage, are assessed in comparison to an existing method [16]. Note, the above mentioned sensor errors and compensation strategies for magnetic disturbance effects are not in the focus here. In the following, Section 1.1 reviews state-of-the-art sensor fusion methods, which provide the basis for the proposed set of methods. Section 1.2 shortly summarizes state-of-the-art methods for obtaining the addressed calibration parameters, which provides an indication of typical error ranges. The contributions of this work are then detailed in Section 1.3.

1.1. Sensor Fusion Methods

Based on the assessed literature, sensor fusion methods are in this work mainly distinguished through the chosen biomechanical model representation and the method for solving the resulting estimation problem. They are also categorized w.r.t. magnetometer usage.

A widespread and efficient representation of a kinematic body model is via a kinematic chain, e.g., [13,16,18,20]. Here, the global orientation and position of the root segment, as well as the relative orientations between segments (cf. Figure 1a), i.e., the joint angles, are modeled as estimation variables. Orientations are typically parametrized in a minimal way, e.g., through Euler angles [13] or Denavit-Hartenberg (DH) coordinates [16,20]. Kinematic chain models offer several advantages: They can be used for predicting body accelerations at the IMU positions, which aids the separation from accelerations due to gravity. In case of a minimal rotation parametrization and, if the segment coordinate systems are sufficiently well aligned with the anatomical rotation axes, restricted joint DoFs can be easily modeled by omitting single angular DoFs. Moreover, relative joint angles are the variables of interest for most applications, e.g., [12,13,14]. At the same time, minimal orientation parametrizations suffer from singularities [26].

Another approach found in literature chooses a redundant parametrization for the biomechanical model by representing each segment with an orientation (mostly via a singularity-free unit quaternion [26]) and position (or velocity [27]) w.r.t. a global frame [17,19,28]. This is subsequently referred to as free segments model (cf. Figure 1b). Conditions from the biomechanical model (e.g., connected segments, restricted rotational DoFs) are then incorporated into the estimation as stochastic constraints, e.g., via virtual measurements [17,19]. In [19] this concept is also extended to the IMU placement. Here, not only global segment kinematics but also the global poses of all IMUs are estimated. The rigid connections between the IMUs and segments, i.e., the I2S calibrations, are then again incorporated via stochastic constraints. Obviously, the definition of a redundant system with stochastic constraints increases the dimensionality of the estimation problem, which results in a higher computational complexity compared to the kinematic chain parametrization (cf.

Table 1. Characteristics of the different sensor fusion methods: n denotes the number of segments and w is the window size used by the Optitracker. All tuning parameters are given in Appendix D. Note, ${\theta }_i$ represent the joint angles estimated by the Chaintracker and ${\Sigma }^{\ddot{\theta }}$ refers to the process noise covariances used in the dynamic model [16].

	Chaintracker (cf. [ 16])	Quattracker segment	Quattracker IMU	Optitracker
Estimation method	EKF	EKF	EKF	WLS
State	${\left({\{{\theta }_i,{\dot{\theta }}_i,{\ddot{\theta }}_i\}}_{i=0}^{n-1}\right)}_t^T$	${\left({\{S_{i,t}^G,{\dot{S}}_{i,t}^G{\ddot{S}}_{i,t}^G,q_{i,t}^{GS},{\omega }_{S,i,t}^{GS},{\dot{\omega }}_{S,i,t}^{GS}\}}_{i=0}^{n-1}\right)}^T$	${\left({\{I_{i,t}^G,{\dot{I}}_{i,t}^G{\ddot{I}}_{i,t}^G,q_{i,t}^{GI},{\omega }_{I,i,t}^{GI}\}}_{i=0}^{n-1}\right)}^T$	${\left({\{{\{S_{i,t}^G,q_{i,t}^{GS},I_{i,t}^G,{\dot{I}}_{i,t}^G,q_{i,t}^{GI}\}}_{i=0}^{n-1}\}}_{t=0}^{w-1}\right)}^T$
Dimensions (state s, meas. vector k)	$x\in {\mathbb{R}}^{9n}$, $s=9n$, $k=7$	$x\in {\mathbb{R}}^{19n}$, $s=19n$, $k=7$	$x\in {\mathbb{R}}^{16n}$, $s=16n$, $k=7$	$x\in {\mathbb{R}}^{17n\times w}$, $s=17nw$
Motion model	1D const angular acc	3D const angular & linear acc	3D const angular vel; 3D const linear accel	IMU control input
Tuning parameters	${\Sigma }^{\ddot{\widehat{\theta }}},{\Sigma }^a,{\Sigma }^{\omega },{\Sigma }^m$	${\Sigma }^{\dot{\widehat{\omega }}},{\Sigma }^{\widehat{p}},{\Sigma }^a,{\Sigma }^{\omega },{\Sigma }^m,{\Sigma }^p,{\Sigma }^G$	${\Sigma }^{\widehat{\omega }},{\Sigma }^{\widehat{p}},{\Sigma }^a,{\Sigma }^{\omega },{\Sigma }^m,{\Sigma }^p,{\Sigma }^G$	${\Sigma }^{cq}$, ${\Sigma }^{cp}$, ${\Sigma }^{\widehat{p}}$, ${\Sigma }^{\dot{\widehat{p}}}$, ${\Sigma }^{\widehat{q}}$, ${\Sigma }^G,{\Sigma }^m,{\Sigma }^{q0},{\Sigma }^p$
Complexity	$\mathcal{O}(k^{2.4}+s^2)$ [49]	$\mathcal{O}(k^{2.4}+s^2)$	$\mathcal{O}(k^{2.4}+s^2)$	$\mathcal{O}(s^3)$ (Gauss Newton method)
Biomech. model	chain	free segments	free segments	free segments
State coordinate system	segment centered	segment centered	IMU centered	IMU and segment centered

). At the same time stochastic constraints consider errors in the associated biomechanical model assumptions and parameters.

Along with the biomechanical model representation there are also different methods for solving the resulting estimation problem. Focusing on online applications the most widespread methods are based on recursive filters, e.g., on an extended Kalman filter (EKF) [9,16,27]. The latter approximates nonlinear motion and measurement equations using a first-order Taylor expansion around the current estimate [29]. The unscented Kalman filter (UKF) [6,20] uses sigma-point approximations [30] that are related to the second-order moments [31]. Note that the UKF does not necessarily yield better performance than the EKF, which has been investigated for the inertial tracking of an arm in [20]. In contrast to these filtering approaches, the work of [19] proposes an optimization-based method to inertial body motion tracking. Here, an offline maximum a posteriori smoothing estimate of the segment kinematics and IMU poses is obtained from a sequence of inertial measurements by solving a global constrained weighted least squares (WLS) problem using an infeasible start Gauss-Newton method. In [32], the estimation problem is decomposed into small subproblems exploiting its sparsity structure. This leads to a distributed, but still offline version of the method, which allows for making use of multiple processors. Using an optimization formulation better accounts for nonlinearities through repeated linearizations. Moreover, it enables the incorporation of global constraints and non-Gaussian noises. Note, the EKF can be considered as a special case of one optimization step [33]. Obviously, an optimization-based estimate is, compared to a filtering approach, computationally more expensive (cf. Table 1). Moreover, the methods in [19,32] are offline, since all measurements need to be collected before processing.

As mentioned above, inertial sensors are typically combined with magnetometers in order to compensate for heading drift resulting from the integration of noisy and biased gyroscope measurements. To this end, the work of [10] provides a survey of different techniques for dealing with magnetic disturbances. From the above referenced methods [18,19,20,27,32] work magnetometer-free and can thus be considered independent of a static magnetic field.

1.2. Calibration Methods

In order to derive I2S orientations, calibration procedures, which require the user to precisely perform predefined static poses [16,17,34] or movements [35,36], are typically used. In a recent study [22] different established calibration methods were validated against an optical reference system based on ten healthy subjects instructed by three operators. The study reports trueness (root mean squared error w.r.t. the optical reference) in the range ${[8;26]}^{\circ }$ and precision (reproducibility) in the range ${[5;10]}^{\circ }$, which attests that potentially large (human-induced) errors w.r.t. the I2S orientations have to be expected, even when subjects are instructed when performing the calibration. Concerning the I2S positions, there are no established calibration methods. In [37], an offline least squares estimator is proposed, which obtains the position of a ball-and-socket joint in the reference frame of two IMUs attached to the adjacent segments from arbitrary motion. The calibration method is used in [18] as basis for knee flexion/extension and ankle plantar/dorsiflexion angle estimation. On simulated data from a three segment model, the precision of the position calibration method is reported with less than three percent difference to the true values. On real data from one subject, the repeatability is reported with variations by about $\pm 0.01$ m. Manually measuring the I2S positions or making assumptions concerning their relative positions w.r.t. the associated segment seem to be the most widespread methods [16,17,19,27,38]. Segment lengths are often also manually measured, based on anatomical landmarks, or they are derived from a measured or assumed height, based on anthropometric tables [16,17,19]. Both methods likely lead to errors in the order of several centimeters, given the difficulty in precisely locating the joint centers, even based on anatomical landmarks [39].

Note, besides the above mentioned position self-calibration method [37], Taetz et al. [40] propose an online capable system, an extension of the optimization-based approach presented in this work, where both I2S positions and orientations are estimated simultaneously with the segment poses. On simulated data from a two segment model, an average precision in the order of sub-degrees for the estimated I2S orientations and in the order of $0.01$ m for the estimated I2S positions is reported. On real data from the lower body of one subject, the repeatability of the I2S orientation estimation is noted with a variation of below $2^{\circ }$ for four IMUs mounted on the upper and lower legs. The estimated I2S positions seemed to be more accurate for the upper legs (repeatability < 0.01 m) than for the lower legs (<0.06 m), which is explained by the low amount of motion variability in the lower legs. Hence, while self-calibration methods appear as a promising way for reducing human-induced calibration errors, they are subject of current research and seem to be not yet in the state of being widely used.

1.3. Contributions

Given that calibration methods based on static poses, movements and manual measurements are still the most widely used, potentially large calibration errors have to be expected, in particular when calibration is performed by inexperienced users. A practical inertial body tracking method should therefore robustly handle such errors. Showing reduced dependence on magnetometer usage also represents an important aspect of robustness in man made environments. See, for instance, [41] for an evaluation of indoor magnetic distortion effects on gait analysis performed with wearable inertial sensors. In this respect, the major contributions of this work are:

• The development of two EKF-based methods with different state-space models, which use the free-segments representation, inspired by [17]. These are subsequently denoted Quattracker IMU and Quattracker segment. Here, rotations are represented through unit quaternions.
• The development of an online capable version of the optimization-based method in [19], based on sliding window optimization. The method is subsequently denoted Optitracker.
• A performance comparison of the new/adapted methods with an existing EKF-based method that uses the kinematic chain representation and DH coordinates to represent joint angles [16]. Performance is measured in terms of angular error statistics on complex (i.e., simultaneous variations in all joint DoFs) moderate and fast human motion (real and simulated) and on artificially simulated complex motion from a case study. In particular, the influence of the selected model calibration errors, i.e., I2S calibration and segment length errors, on the performances of the different methods and their dependence on magnetometer usage are assessed.

Based on the results from these studies, useful considerations concerning the usage of the different tested methods are provided. In the following, the notation, methods and evaluation setup are introduced in Section 2. The results are summarized in Section 3. Section 4 discusses these results and draws conclusions.

2. Materials and Methods

2.1. Notation

Let $A,B$ be two Cartesian coordinate systems (coordinate frames), then $R^{AB}\in SO(3)$ denotes the orientation of frame B in frame A and $A^B\in {\mathbb{R}}^3$ denotes the translation of frame B in frame A. The rotation matrix $R^{AB}$ and its unit quaternion representation $q^{AB}$ are used interchangeable. Moreover, $H^{AB}$ denotes a homogeneous transformation comprising both, orientation and position $\{R^{AB},A^B\}$. To switch between the representations, $rot(H^{AB})R^{AB}$ extracts the rotation matrix from a homogeneous transformation. The translation is extracted with $trans(H^{AB})A^B$. The inverse rotation is denoted ${(q^{AB})}^{*}{q}^{BA}$, ${(R^{AB})}^T=R^{BA}$. The angular velocity of frame B w.r.t. A in frame A is denoted ${\omega }_A^{AB}$. The cross product $v\times w$ with $v,w\in {\mathbb{R}}^3$ can be rewritten as a matrix multiplication $S(v)w$, with $S(v)$ being a skew-symmetric matrix. Let $a_t\in {\mathbb{R}}^m$ be a time dependent variable of dimension m, then ${\dot{a}}_t$ denotes the first and ${\ddot{a}}_t$ the second time derivative. Process and measurement noises are assumed to be additive and Gaussian distributed with zero mean. Measurement noises are denoted $e_t^X\sim \mathcal{N}(0,{\Sigma }^X)$, with ${\Sigma }^X$ being the noise covariance matrix, while process noises are denoted $e_t^{\widehat{X}}\sim \mathcal{N}(0,{\Sigma }^{\widehat{X}})$. Acceleration, angular velocity and magnetic field measurements from IMU i at time t are denoted $y_{i,t}^a,y_{i,t}^{\omega },y_{i,t}^m\in {\mathbb{R}}^3$, respectively.

2.2. Biomechanical Model Representations

The two biomechanical model representations, i.e., the kinematic chain model and the free segments model (see Figure 1), are now formalized using the introduced notation. Both representations share a global position-less coordinate system G, whose x-axis is aligned with the local magnetic north and whose z-axis points up, opposite gravity. Each segment i has a local coordinate frame $S_i$ and is assumed to have one IMU attached to it. Thus, for notational brevity, segment and attached IMU share one index. Each IMU i has a local coordinate frame $I_i$, in which the measurements are represented. Moreover, $S_i$ and $I_i$ are rigidly connected via the I2S orientation $q_i^{SI}$ and position $I_i^S$ (or $H_i^{SI}$ for the chain model). The kinematic chain model, illustrated in Figure 1a, defines a hierarchical order of transformations. In this work, the definitions from [16] are used, which are based on DH transformations [42]. For the sake of completeness, the transformations and equations for the model construction are provided in Appendix A. Further details can be found in [16]. The free segments model is composed of a set of segments $S_i\in \mathbb{S}$, a set of IMUs $I_i\in \mathbb{I}$ and a set of joints $J_k\in \mathbb{J}$. For each intermediate joint $J_k$, the set of connected segments is denoted $S_i,S_j\in {\mathbb{S}}_{J_k}$. Each segment $S_i$ and IMU $I_i$ is represented in G with an orientation and translation {$q_i^{GS}$, $S_i^G$} and {$q_i^{GI}$, $I_i^G$}, respectively. Figure 1b illustrates this model. Note, each segment also has a start and an endpoint defined in the local segment coordinate system ($p_i^S$, $p_j^S$ in the figure), which are used in Equation (10).

2.3. EKF-Based Methods

The Chaintracker, the Quattracker IMU and the Quattracker segment use an EKF for parameter estimation. The former is taken from [16] and is enhanced with an initialization method described in Section 2.3.5. The other methods are inspired by the work of [17,19]. Note, the work of [17] provides only a general concept for combining a free segments model with biomechanical constraints. The formalizations with different state-space models leading to the Quattracker IMU and Quattracker segment were contributed as part of this manuscript.

2.3.1. Measurement Models

For IMU $I_i$, the accelerometer measurement model at time t is:

$$y_{i,t}^a={(R_{i,t}^{GI})}^T\left({\ddot{I}}_i^G-g^G\right)+e_{i,t}^a$$

(1)

Here, ${\ddot{I}}_i^G$ and $g^G$ denote body acceleration and acceleration due to gravity, respectively, both given in the global frame.

The gyroscopes measure the angular velocities w.r.t. the global frame, transformed into the IMU frame, resulting in the following gyroscope measurement model:

$$\begin{array}{ccc}\hfill y_{i,t}^{\omega }& =& {\omega }_{I,i,t}^{GI}+e_{i,t}^{\omega }\hfill \end{array}$$

(2)

Note, gyroscope and accelerometer bias models are not in the focus of this paper, but can be easily added, see e.g., [19,43].

The magnetometer measurement model was chosen to only have an effect on the estimated yaw direction. It omits information concerning the local dip angle. This is a common way of reducing the influence of magnetic disturbances [10]. Let ${\widehat{y}}_{i,t}^m$ be the normalized magnetometer measurement, then the model can be written as:

$$0=\mathrm{atan}2\left[\frac{{(R_{i,t}^{GI}{\widehat{y}}_{i,t}^m)}_y}{{(R_{i,t}^{GI}{\widehat{y}}_{i,t}^m)}_x}\right]+e_{i,t}^m$$

(3)

Here, ${(·)}_x$, ${(·)}_y$ denote the x- and y-component of the vector. Hence, the model considers the angular deviation of the magnetometer measurement, as transformed into the global frame and projected into the horizontal plane, from the global x-axis.

In all models, the three variables ${\ddot{I}}^G$, ${\omega }_I^{GI}$ and $R^{GI}$ need to be extracted from the state space, as defined in the following.

2.3.2. State Spaces

The Quattracker IMU uses a constant angular velocity and a constant acceleration model [43]. Thus, the angular velocity and the linear acceleration are required in the state. Assuming a rigid and known I2S calibration, the IMU positions and orientations are sufficient to reconstruct the corresponding segment pose. Hence, for n segments the state at time t comprises:

$$x_t={\left({\{I_{i,t}^G,{\dot{I}}_{i,t}^G{\ddot{I}}_{i,t}^G,q_{i,t}^{GI},{\omega }_{I,i,t}^{GI}\}}_{i=0}^{n-1}\right)}^T$$

(4)

The Quattracker segment, similarly to the Chaintracker, represents the kinematic variables in the segment frame. This leads to a dynamic model assuming constant linear and angular acceleration, while the angular acceleration is required as part of the state:

$$x_t={\left({\{S_{i,t}^G,{\dot{S}}_{i,t}^G{\ddot{S}}_{i,t}^G,q_{i,t}^{GS},{\omega }_{S,i,t}^{GS},{\dot{\omega }}_{S,i,t}^{GS}\}}_{i=0}^{n-1}\right)}^T$$

(5)

Note, the introduction of $\dot{\omega }$ is a consequence of the need for ${\ddot{I}}^G$ in Equation (1), since:

$$\begin{array}{ccc}\hfill I_{i,t}^G& =& S_{i,t}^G+R_{i,t}^{GS}{I}_{i,t}^S\hfill \end{array}$$

(6a)

$$\begin{array}{ccc}\hfill \Rightarrow {\dot{I}}_{i,t}^G& =& {\dot{S}}_{i,t}^G+{\dot{R}}_{i,t}^{GS}{I}_{i,t}^S+R_{i,t}^{GS}\underset{=0}{\underbrace{{\dot{I}}_{i,t}^S}}\hfill \end{array}$$

(6b)

$$\begin{array}{ccc}\hfill {\dot{I}}_{i,t}^G& =& {\dot{S}}_{i,t}^G+R_{i,t}^{GS}S({\omega }_{S,i,t}^{GS})I_{i,t}^S\hfill \end{array}$$

(6c)

$$\begin{array}{ccc}\hfill \Rightarrow {\ddot{I}}_{i,t}^G& =& {\ddot{S}}_{i,t}^G+{\dot{R}}_{i,t}^{GS}S({\omega }_{S,i,t}^{GS})I_{i,t}^S+R_{i,t}^{GS}S({\dot{\omega }}_{S,i,t}^{GS})I_{i,t}^S\hfill \end{array}$$

(6d)

2.3.3. Dynamic Models

The dynamic model for the Quattracker IMU is:

$$x_{t+T}=\left[\begin{array}{c}{I}_{i,t+T}^G\\ {\dot{I}}_{i,t+T}^G\\ {\ddot{I}}_{i,t+T}^G\\ q_{i,t+T}^{GI}\\ {\omega }_{I,i,t+T}^{GI}\\ ⋮\end{array}\right]=\left[\begin{array}{c}{I}_{i,t}^G+T{\dot{I}}_{i,t}^G+\frac{T^2}{2}{\ddot{I}}_{i,t}^G\\ {\dot{I}}_{i,t}^G+T{\ddot{I}}_{i,t}^G\\ {\ddot{I}}_{i,t}^{GI}+T{e}_{i,t}^{\widehat{p}}\\ q_{i,t}^{GI}\odot exp(\frac{T}{2}{\omega }_{I,i,t}^{GI})\\ {\omega }_{I,i,t}^{GI}+T{e}_{i,t}^{\widehat{\omega }}\\ ⋮\end{array}\right]$$

(7)

where ⊙ is the quaternion product, $exp$ denotes the quaternion exponential and T is the sampling time. Note, the vertical dots indicate that these variables are given for each segment $i\in \{0,\dots ,n-1\}$.

Equation (7) is built similarly for the Quattracker segment, by applying the following two modifications: (1) replace $I^G$ with $S^G$ and (2) consider $\dot{\omega }$ in the rotational update by utilizing Equation (B4) (Appendix B). Thus, Equation (7) can be extended as follows:

$$x_{t+T}=\left[\begin{array}{c}{S}_{i,t+T}^G\\ {\dot{S}}_{i,t+T}^G\\ {\ddot{S}}_{i,t+T}^G\\ q_{i,t+T}^{GS}\\ {\omega }_{S,i,t+T}^{GS}\\ {\dot{\omega }}_{S,i,t+T}^{GS}\\ ⋮\end{array}\right]=\left[\begin{array}{c}{S}_{i,t}^G+T{\dot{S}}_{i,t}^G+\frac{T^2}{2}{\ddot{S}}_{i,t}^G\\ {\dot{S}}_{i,t}^G+T{\ddot{S}}_{i,t}^G\\ {\ddot{S}}_{i,t}^{GI}+T{e}_{i,t}^{\widehat{p}}\\ q_{i,t}^{GS}\odot exp\left(\frac{T}{2}{\omega }_{S,i,t}^{GS}+\frac{T^2}{4}{\dot{\omega }}_{S,i,t}^{GS}\right)\\ {\omega }_{S,i,t}^{GS}+T{\dot{\omega }}_{S,i,t}^{GS}\\ {\dot{\omega }}_{S,i,t}^{GS}+e_{i,t}^{\dot{\widehat{\omega }}}\\ ⋮\end{array}\right]$$

(8)

The following equations obtain the quantities required for the measurement models (as defined in Section 2.3.1) from the above state using the I2S calibrations:

$$\begin{array}{cc}\hfill {\ddot{I}}_{i,t}^G& ={\ddot{S}}_{i,t}^G+R_{i,t}^{GS}\left[S({\omega }_{S,i,t}^{GS})S({\omega }_{S,i,t}^{GS})+S({\dot{\omega }}_{S,i,t}^{GS})\right]I_i^S\hfill \end{array}$$

(9a)

$$\begin{array}{cc}\hfill q_{i,t}^{GI}& =q_{i,t}^{GS}\odot q_i^{SI}\hfill \end{array}$$

(9b)

$$\begin{array}{cc}\hfill {\omega }_I^{GI}& =q_i^{IS}\odot {\omega }_{S,i,t}^{GS}\odot q_i^{SI}={(R_{i,t}^{SI})}^T{\omega }_{S,i,t}^{GS}.\hfill \end{array}$$

(9c)

2.3.4. Constraints

Constraints within an EKF are commonly implemented as measurement models [44,45]. In this work, in order to reduce drifts, constraints are used to ensure the segments staying connected at the joints. The following models are reformulated biomechanical constraints from [19]. Let $S_{i,t}^G,S_{j,t}^G$ be the global position of two connected segments $S_i,S_j\in {\mathbb{S}}_{J_k}$ at time t. Let $R_{i,t}^{GS}$ and $R_{j,t}^{GS}$ be their orientations. Given $p_i^S\in {\mathbb{R}}^3$, an endpoint of $S_i$, and $p_j^S\in {\mathbb{R}}^3$, the corresponding point in $S_j$ (cf. Figure 1b), the measurement model that ensures connection of the two segments is:

$$0=S_{i,t}^G+R_{i,t}^{GS}{p}_i^S-(S_{j,t}^G+R_{j,t}^{GS}{p}_j^S)-e_{i,t}^p$$

(10)

By exchanging the point of segment j by a constant point in the global system $P^G$, the model:

$$0=S_{i,t}^G+R_{i,t}^{GS}{p}_i^S-P^G-e_{i,t}^G,$$

(11)

ensures that the point $p_i^S$ stays close to the global point $P^G$ for all time steps. This is used in the evaluation for fixing the position of the root segment’s origin.

2.3.5. Initialization

Given the first set of measurements, the TRIAD method [46] is used to determine the initial IMU orientations $q_{i,0}^{GI}$. For the Quattracker IMU the obtained orientations can be introduced into the state directly, whereas the Quattracker segment requires these to be transformed into the segment orientations using the I2S calibrations. After initializing the orientations, in a second step, it is ensured, that Equation (10) is fulfilled for the correctly rotated segments. The initial angular and linear velocities and accelerations are set to zero. The state of the Chaintracker comprises the variable angles and angle derivatives of the given kinematic chain (Appendix A). The angle derivatives are initialized with zero, while the angles are calculated from the initial quaternions $q_{i,0}^{GI}$ via inverse kinematics. This is done by minimizing:

$$\underset{x_c}{argmin}\sum _{i=0}^{n-1}{\parallel (R_{i,0}^{GI}rot{(f_i(x_c))}^T-I_{3\times 3})\parallel }_F$$

(12)

where n is the number of IMUs, ${\left|\right|·\left|\right|}_F$ denotes the Frobenius norm and $rot{(f(x_c))}^T$ extracts the global orientations from the chain state $x_c$ (Appendix A).

2.4. Sliding Window Optimization

While the presented EKF approaches only maintain a state for a single time step and process one IMU data set at a time, the optimization-based approach uses windows or batches of IMU data sets to compute a corresponding batch of states (see Equation (13)). The method is based on [19], however, extends this offline optimization approach to an online-capable method by introducing the window mechanism and appropriately adapting the cost function (cf. Section 2.4.2). As suggested in [19], both the IMU and segment poses are estimated in the state. This relaxes the rigidity assumption concerning the I2S calibrations. Moreover, in contrast to the EKF approaches, the dynamic model takes the IMU data as control input (see Appendix C). This results in the IMU velocities being estimated in the state, while the need for estimating the linear accelerations and angular velocities is avoided. The state is composed as follows:

$$x^b={\left({\{{\{I_{i,t}^G,{\dot{I}}_{i,t}^G,q_{i,t}^{GI},S_{i,t}^G,q_{i,t}^{GS}\}}_{i=0}^{n-1}\}}_{t=0}^{w-1}\right)}^T$$

(13)

where w is the window size, n the number of segments and b is the batch number. For every optimization step, an overlap $ov$ with $w>ov\ge 1$ is defined. Setting $ov=1$ results in a delayed tracking, where the delay depends on the setting of w, while $ov=w-1$ leads to online tracking using moving horizon estimation. Here, for each new time step, the latest IMU data is added and the oldest is discarded.

The optimization-based approach comes down to solving a weighted least-squares problem for each batch of data, where the different terms in the cost function are obtained by rearranging the stochastic model equations so that the noises are isolated on the left-hand side:

$$\begin{array}{c}\hfill \underset{x^b}{min}\sum _{t=0}^{w-1}\sum _{i=0}^{n-1}\left(\underset{12\mathrm{S}\text{calibrations},\text{Appendix C}}{\underbrace{{∥e_{i,t}^{cq}∥}_{{({\Sigma }^{cq})}^{-1}}^2+{∥e_{i,t}^{cp}∥}_{{({\Sigma }^{cp})}^{-1}}^2}}\right.+\underset{\text{motion}\text{model},\text{Appendix C}}{\underbrace{{∥e_{i,t}^{\widehat{p}}∥}_{{({\Sigma }^{\widehat{p}})}^{-1}}^2+{∥e_{i,t}^{\dot{\widehat{p}}}∥}_{{({\Sigma }^{\dot{\widehat{p}}})}^{-1}}^2+{∥e_{i,t}^{\widehat{q}}∥}_{{({\Sigma }^{\widehat{q}})}^{-1}}^2}}\\ \hfill \left.+\underset{\text{prior},\text{Equation}(11)}{\underbrace{\parallel e_{i,t}^G{\parallel }_{{({\Sigma }^G)}^{-1}}^2}}+\underset{\text{magnetometer}\text{model},\text{Equation}(3)}{\underbrace{\parallel e_{i,t}^m{\parallel }_{{({\Sigma }^m)}^{-1}}^2}}\right)\\ \hfill +\sum _{i=0}^{n-1}\left(\underset{\text{initialization},\text{Section}\mathrm{2.4.2}}{\underbrace{\parallel e_{i,t}^{q0}{\parallel }_{{({\Sigma }^{q0})}^{-1}}^2}}\right)+\sum _{t=0}^{w-1}\sum _{J_k\in \mathbb{J}}\left(\underset{\text{biomechanical}\text{constraints},\text{Section}\mathrm{2.4.1}}{\underbrace{\parallel (e_{k,t}^p){\parallel }_{{({\Sigma }^p)}^{-1}}^2}}\right)\end{array}$$

(14)

The following sections explain the individual terms in more detail. For the sake of completeness, the terms already proposed in [19] are provided in the Appendix C. Note, compared to [19], the prior and magnetometer model term were added in order to enable fair comparison with the EKF-based methods.

2.4.1. Biomechanical Constraints

The biomechanical constraints (connected joints) are based on the re-arrangement of Equation (10). In contrast to [19], where these are integrated as hard constraints, these are here included as terms into the cost function. Thus, an unconstrained weighted least squares problem is obtained, which can be solved using standard nonlinear least-squares methods, for instance the Gauss Newton or Levenberg Marquardt method [47,48].

2.4.2. Initialization

The initialization of the IMU orientations can be obtained by minimizing $\forall I_i\in \mathbb{I}$

$$e_{i,t}^{q0}=\left\{\begin{array}{cc}2log\left(q_{i,w-1}^{IG,b-1}\odot q_{i,0}^{GI,b}\right)\hfill & \text{for}b>0\hfill \\ 2log\left(q_{i,init}^{IG}\odot q_{i,0}^{GI,b}\right)\hfill & \text{else}\hfill \end{array}\right.$$

(15)

This penalizes sudden changes of the estimated IMU orientations for the overlap of $b-1$ and b, which is a required adaptation for the proposed online-capable method compared to [19]. Note, in a moving horizon context this term corresponds to the arrival cost for the variables. Here, $log$ denotes the quaternion logarithm. Note, for $b=0$, similarly to the EKF methods, the initial quaternions $q_{i,init}^{GI}$ are obtained using the TRIAD algorithm.

2.5. Summary and Overview

After the different methods have been described in detail, Table 1 provides a comparative summary of their characteristics.

It also adds information concerning their computational complexities, which mainly depend on the number of estimation variables per time step and the solution method. Also note that a higher redundancy in the formulation of the estimation problem results in an increase in the number of tuning parameters.

The Chaintracker, in contrast to all other methods, only models rotational quantities in the state. This is due to using the kinematic chain model where positions and translations are implicitly given. All EKF-based methods assume the I2S calibrations being error-free, while the Optitracker assumes Gaussian distributed zero-mean errors being present in the I2S calibrations. The Quattracker IMU and the Quattracker segment differ in the coordinate frame, in which the state variables are given. The former is IMU centered, while the latter, similarly to the Chaintracker, is segment centered. The segment centered formulation results in the need for keeping the segments’ angular accelerations in the state and, thus, taking these into account in the dynamic model. Hence, the Quattracker segment and the Chaintracker only differ in the biomechanical model representation (kinematic chain vs. free segments) and the orientation parametrization (quaternions vs. joint angles). The Optitracker, in contrast to the EKF-based methods, keeps both the segment poses and the IMU poses in the state. Moreover, it takes the IMU data as control input to the dynamic model. With its comparably high level of redundancy in both spatial and temporal dimension, as well as, its iterative estimation method, the Optitracker has a significantly higher computational complexity than the other methods.

2.6. Evaluation Setup

The overall goal of the evaluation was to answer the question, which of the proposed sensor fusion schemes (Quattracker segment, Quattracker IMU, Optitracker), in comparison to the existing method (Chaintracker, here considered as baseline), provides the best basis for developing a truly robust online inertial body tracking method. Based on the argumentation in Section 1, sensitivity to the selected model calibration errors and dependence on magnetometer usage was in the focus. The evaluation was performed as a case study with a kinematic model comprising three rigid segments with known lengths, which are connected through two three DoF joints. The root joint was assumed fixed in space, with three rotational DoFs. In terms of motion sequences, moderate and fast motion was considered, in order to compare the performances of the different sensor fusion methods w.r.t. the above mentioned criteria also in relation to the motion agility (cf. [20]). General complex (i.e., containing simultaneous variations in all joint DoFs) motion was chosen instead of planar motions or specific activities, since the focus was on the comparison of the different sensor fusion methods under general conditions, rather than on their ability to represent specific motion patterns. Further explanations are provided below. Within this case study, we investigated two scenarios: (1) a real data scenario and (2) a simulation scenario with systematically introduced calibration errors. Both were investigated under two test configurations: (1) using magnetometer information (w/mag) and (2) dropping the magnetometer information, i.e., dropping the magnetometer measurement model in Equation (18) for all sensor fusion methods (w/o mag). The two scenarios, together with the used data sequences, as well as, the measures of performance are described in the following.

2.6.1. Real Data Scenario

The goal of the real data scenario was to assess and compare the overall performances of the different sensor fusion methods under realistic conditions. Here, different error sources are unavoidably present and cannot be reliably separated, as also discussed in [9]. The major ones are therefore characterized in Section 2.6.2. Hence, the real data setup served for gathering tendencies concerning the overall robustness of the different sensor fusion methods under realistic conditions, while the influence of the considered model calibration errors were assessed in simulation (cf. Section 2.6.3). The results of the real data scenario are summarized in Section 3.1.

For data collection, the following systems were used: the NaturalPoint OptiTrack system with 12 Prime 13 cameras, operated with the Motive software (Version 1.8) (NaturalPoint, Inc., Corvallis, OR, USA) [50] and the XSens Link system, operated with the MVN Studio BIOMECH software (Version 4.1) (Xsens, Enschede, The Netherlands) [51]. The former served as reference, while the latter provided the real IMU data. Both systems recorded data at 120 Hz.

A healthy male (age: 30 years, height: $1.76$ m) was equipped with three IMUs at the right upper arm, forearm and hand (cf. Figure 2). Each IMU was mounted into a special casing with optical markers (4 mm diameter) attached to it, in order to ensure a rigid connection between IMU and marker rigid body. Additional markers (12 mm diameter) were placed at anatomical landmarks in order to enable skeleton fitting using the Motive software (upper body model, 25 markers). The cameras were arranged in a small volume in order to enable continuous tracking of the small markers.

With this setup, two evaluation data sequences with complex moderate and fast human motion were recorded. For this, the subject was asked to perform a movement, which introduces simultaneous variations in all DoFs of the right arm, while keeping the shoulder stationary, once at moderate speed and once at fast speed with fast speed changes. The movements also had to be performed anterior to the frontal plane and with the hand roughly below the shoulder height in order to be well captured by the optical system. The resulting data sequences, subsequently denoted real-slow and real-fast, both contain eight-shaped movements at respective speeds, smooth parts reminding of reaching and steering in the case of real-slow, and agile parts with quick starts and stops, as well as, parts reminding of sportive movements, such as boxing, in the case of real-fast. In Figure 3, the two data sequences are illustrated in terms of Euler angles per time step and respective ranges of motion for each joint DoF. The ranges of motion are comparable between the sequences, besides sim-fast covering smaller ranges for rotations around the y-axes of shoulder and elbow (external rotations of upper and forearm). The ranges of motion for the hand DoFs are the smallest among the considered joints. Videos showing re-simulations of real-slow and real-fast (called sim-slow and sim-fast) are available as supplementary files.

For each time step t, the data sequences contain the orientations and positions of the marker rigid bodies (as generated with the Motive software) in the world reference frame O of the optical system, the joint positions and orientations of the fitted skeleton and the IMU data, ${\{y_{i,t}^{\omega },y_{i,t}^a,y_{i,t}^m\}}_{i=0}^2$. The joint data was only used to obtain the I2S calibrations (see below). In addition to the Euler angles, the measured peak accelerometer and angular velocity 2-norms are shown in

Table 2. Measured/simulated instantaneous peak acceleration (Acc) and angular velocity (Gyr) 2-norms for real-slow, real-fast, sim-slow and sim-fast.

Mode	Sequence →	Slow		Fast
	Sensor →	Acc(m/s${}^2$)	Gyr (${}^{\circ }/\mathrm{s}$)	Acc (m/s${}^2$)	Gyr (${}^{\circ }/\mathrm{s}$)
Real	$I_0$	$13.44$	$178.68$	$26.88$	$457.05$
	$I_1$	$16.87$	$394.68$	$75.24$	$957.29$
	$I_2$	$18.74$	$468.55$	$95.45$	$1031.31$
Re-simulated	$I_0$	$18.05$	$187.07$	$25.70$	$429.75$
	$I_1$	$22.07$	$395.27$	$77.20$	$968.35$
	$I_2$	$22.86$	$440.23$	$100.55$	$1030.76$

.

Similarly to [52], temporal synchronization of optical and inertial data was done by maximizing the correlation between the angular velocities measured by the IMUs and the angular velocities derived from the marker rigid body orientations. The latter were calculated according to [53].

For spatial synchronization of the orientations obtained from the marker rigid bodies, $q_{i,t}^{OB}$, and the segment orientations estimated by the different sensor fusion methods, $q_{i,t}^{GS}$, the marker rigid bodies were transformed into the IMU reference frames using:

$$\begin{array}{c}\hfill {\overline{q}}_{i,t}^{GS}=q^{GO}\odot q_{i,t}^{OB}\odot q_i^{BI}\odot q_i^{IS}\end{array}$$

(16)

Here, $q^{GO}$ and $q_i^{BI}$ denote the hand-eye orientations, which are the relative orientations aligning the global optical frame O with the global inertial frame G, and the marker rigid body frames $B_i$ with the IMU frames $I_i$. These were estimated from a data sequence with smooth motion, where the three IMUs were mounted on a stick, using the approach of [54]. Since the alignment approach in [54] is based on rotation sequences, the fused IMU orientations provided by the Link IMUs, $q_{i,t}^{GI}$, rather than the measured IMU data, and the orientations of the marker rigid bodies, $q_{i,t}^{OB}$, were used for calculation. Moreover, $q_i^{IS}$ refers to the I2S orientations. The I2S calibrations were calculated from a static recording of the subject being in a T-pose. For this, the fitted skeleton (i.e., the fitted 3D joint positions and orientations), the marker rigid body poses and the hand-eye orientations were used. The segment lengths were also derived from the fitted skeleton with: $\parallel S_0\parallel =0.28$ m, $\parallel S_1\parallel =0.23$ m, $\parallel S_2\parallel =0.08$ m (including only the palm). The resulting biomechanical model is shown in Figure 2 (left).

2.6.2. Real Data Scenario: Discussion of Major Error Sources

This section provides a characterization of the major error sources in both the optical reference data and the IMU data. Such errors are unavoidably present and have to be taken into account for the analysis provided in Section 3.1. Equation (16) shows the different components involved into obtaining the reference segment orientations, which are subsequently used for measuring the performances of the different methods (cf. Section 2.6.4). Here, the errors in the orientations of the marker rigid bodies, $q_{i,t}^{OB}$, as obtained from the optical system, cannot be estimated. Note, the point tracking is claimed to provide sub-millimeter precision. The angular errors contributed through the hand-eye orientations, $q^{GO}$, $q_i^{BI}$, can be quantified for each IMU i to some extend based on the residual errors as calculated from each time step t of the data sequence used for calibration:

$$E_{i,t}^{he}=2\mathrm{acos}(q^{GO}\odot q_{i,t}^{OB}\odot q_i^{BI}\odot q_{i,t}^{IG})$$

(17)

Here, as mentioned above, $q_{i,t}^{IG}$ refers to the orientations obtained from the IMUs. The results are given in

Table 3. Mean (std,max) angular residual errors (cf. Equation (17)) for the hand-eye calibrations of each inertial measurement unit (IMU) $I_i$ as calculated on the data sequence used for calibration.

IMU	Residual Error
$I_0$	$1.{14}^{\circ }(0.{57}^{\circ },3.{89}^{\circ })$
$I_1$	$2.{28}^{\circ }(0.{90}^{\circ },6.{10}^{\circ })$
$I_2$	$1.{51}^{\circ }(0.{77}^{\circ },5.{07}^{\circ })$

.

Concerning the I2S calibrations, $q_i^{IS}$, there are two aspects to consider. First, the errors in the calculated quantities cannot be estimated, since these were derived based on data provided by the optical system (see above). Second, due to soft tissue artifacts, these calibrations are in reality not rigid. In our study, the influence of soft tissue artifacts was not in focus and was minimized by well suited sensor placement, tight fixation and by ensuring a rigid connection between the marker rigid bodies used for calculating the reference segment orientations and the IMUs. Since the proposed sensor fusion methods all assume a kinematic model with rigid I2S calibrations (the Optitracker being the only method, which models Gaussian distributed I2S calibration errors), soft tissue artifacts, which can partly be interpreted as time-dependent changes of the I2S calibrations, still influence the estimated segment orientations, whereas the resulting error cannot be assessed.

The Link IMUs are intrinsically calibrated by the manufacturer, however, small gyroscope biases remain. These were approximated as mean gyroscope values from a static IMU data sequence and were subtracted from the angular velocities of the real-slow and real-fast sequences. The mean 2-norms of the accelerometer measurements under static conditions were $9.74$ m/s${}^2$, $9.76$ m/s${}^2$, $9.77$ m/s${}^2$ for the three IMUs. These deviations from the local gravity strength [55] of $9.81$ m/s${}^2$ were, however, not corrected.

Though no ferromagnetic objects where present in the capturing volume during recording, the measured magnetic field vectors indicate an inhomogeneous magnetic field.

Table 4. Mean (std) of magnetic field vector 2-norms (upper values) and global angular deviations (lower values) for each IMU as calculated from the real data sequences.

	$I_0$	$I_1$	$I_2$
real-slow	$0.92(0.00)$	$0.90(0.01)$	$0.90(0.02)$
	$1.{79}^{\circ }(1.{18}^{\circ })$	$3.{66}^{\circ }(2.{41}^{\circ })$	$4.{36}^{\circ }(2.{89}^{\circ })$
real-fast	$0.92(0.00)$	$0.91(0.01)$	$0.91(0.02)$
	$2.{63}^{\circ }(2.{38}^{\circ })$	$4.{45}^{\circ }(5.{60}^{\circ })$	$4.{59}^{\circ }(6.{04}^{\circ })$

shows the statistics of the magnetic field vector 2-norms and the angular deviations, $\measuredangle (y_{i,t}^m,y_{mean,t})$, in the global frame for real-slow and real-fast. For each time step t, the latter were calculated as the normalized magnetic field measurements transformed into the global frame G using the optical reference data. More specifically, for n IMUs:

$$\begin{array}{ccc}\hfill y_{mean,t}& =& \frac{1}{n}\sum _{i=0}^n{R}_{i,t}^{GI}\frac{y_{i,t}^m}{\parallel y_{i,t}^m\parallel }\hfill \end{array}$$

(18a)

$$\begin{array}{ccc}\hfill \measuredangle (y_{i,t}^m,y_{mean,t})& =& \mathrm{acos}\left(y_{mean,t}{R}_{i,t}^{GI}\frac{y_{i,t}^m}{\parallel y_{i,t}^m\parallel }\right)\hfill \end{array}$$

(18b)

where $R_{i,t}^{GI}$ were extracted from the optical reference data (cf. Equation (16)). Note that the resulting errors include the hand-eye calibration errors (cf. Table 3).

The table indicates stronger magnetic disturbances for $I_1$ and $I_2$ than for $I_0$, which is due to the varying distances of those IMUs from the floor, the latter causing an inhomogeneous magnetic field (as also observed in [56,57]). Since this work investigates the dependence of the proposed sensor fusion methods on magnetometer usage, rather than active compensation of magnetic disturbances, for evaluation, we also tested the different methods on real inertial data in combination with simulated magnetic field vectors (as also done in [9]). The latter were calculated by rotating the x-axis of the global frame G into the frame of each IMU i using:

$$\begin{array}{c}\hfill y_{art,i,t}^m=R_{i,t}^{IG}{(1,0,0)}^T\end{array}$$

(19)

Here, again, the rotations $R_{i,t}^{GI}$ were extracted from the optical reference data (cf. Equation (16)).

2.6.3. Simulation Scenario with Systematically Introduced Model Calibration Errors

The goal of the simulation scenario was to assess the influence of model calibration errors, i.e., I2S calibration and segment length errors, on the different sensor fusion methods, in the absence of other error sources. For this, three simulated data sequences were considered: a simulated data sequence based on an artificially generated smooth, but fast motion and re-simulations of real-slow and real-fast, which both contain human motion, as further explained below. These data sequences are subsequently referred to as sim-fast-artificial, sim-slow and sim-fast, respectively. The results of the simulation scenario are summarized in Section 3.2 and Section 3.3. The sim-fast-artificial sequence was simulated from an arm-like three segment kinematic chain model, which was parametrized using DH coordinates (cf.

Table 5. Denavit-Hartenberg (DH) coordinates for the three segment kinematic chain model used for simulating the sim-fast-artificial data sequence. The angles ${\alpha }_{[0:5]}(t)$ and ${\theta }_{[0:2]}(t)$ are the Degrees of Freedom (DoFs) that are controlled via Equation (20). The Inertial Measurement Unit (IMU)-to-Segment (I2S) positions are given by translations along the bone, in z-direction relative to the segment origins (i.e., DH$(z,0,0,0)$). The initial chain configuration (pointing up opposite gravity) is illustrated on the right.

Segment ($S_i$)	d	a	${\alpha }_j(t)$	${\theta }_j(t)$	IMU	Image
$S_0$	0	0	${\alpha }_0(t)$	$-\frac{pi}{2}$	None
	0	0	${\alpha }_1(t)$	$\frac{pi}{2}$	None
	0	0	0	${\theta }_0(t)$	None
	$0.4$	0	0	0	$z=0.3$
$S_1$	0	0	${\alpha }_2(t)$	$-\frac{pi}{2}$	None
	0	0	${\alpha }_3(t)$	$\frac{pi}{2}$	None
	0	0	0	${\theta }_1(t)$	None
	$0.4$	0	0	0	$z=0.3$
$S_2$	0	0	${\alpha }_4(t)$	$-\frac{pi}{2}$	None
	0	0	${\alpha }_5(t)$	$\frac{pi}{2}$	None
	0	0	0	${\theta }_2(t)$	None
	$0.2$	0	0	0	$z=0.1$

). The segment lengths were $\parallel S_0\parallel =0.4,\parallel S_1\parallel =0.4,\parallel S_2\parallel =0.2$, the I2S positions were $I_0^S={(0,0,0.3)}^T$, $I_1^S={(0,0,0.3)}^T$, $I_2^S={(0,0,0.1)}^T$ (all in m) and the I2S orientations were assumed identity.

In order to animate this kinematic chain, $S_0$ was kept stationary at $W^G={(0,0,0.5)}^T$ and an angle sequence ${\{\varphi \}}_{t=0}^{628}$ was generated using:

$$\varphi (\beta (t))=sin\left(\frac{\beta (t)}{2}\right)sin(\beta (t))\pi$$

(20)

with $\beta (t)=\frac{2\pi }{629}t$. This results in a sampling of one $2\pi$ period with a discretization of roughly $0.01$ radians, where each step is assumed to correspond to a sampling time of $0.01$s. This sequence was used for each rotational DoF ${\alpha }_0,\dots ,{\theta }_2$ and the resulting motion provides the ground truth segment orientations and positions for the sim-fast-artificial sequence. The angle sequence is visualized in Figure 4 and the motion resulting from applying this angle sequence to each rotational DoF is provided as a supplementary video. Based on the segment kinematics, the IMU trajectories are calculated by applying the assumed I2S calibrations. From this, IMU data was obtained using standard data differentiation [26] and, if applicable, by applying realistic sensor noises as detailed in Appendix D. The peak acceleration and angular velocity 2-norms resulting from the above described I2S calibrations for sim-fast-artificial are given in

Table 6. Peak acceleration (Acc) and angular velocity (Gyr) 2-norms for sim-fast-artificial. For all sensors, the values vary smoothly between 0 and the peak values shown in the table.

Sequence →	sim-fast-artificial
Sensor →	Acc (m/s${}^2$)	Gyr (${}^{\circ }/\mathbf{s}$)
$I_0$	$14.03$	$356.90$
$I_1$	$38.16$	$705.20$
$I_2$	$61.90$	$1047.99$

.

The sim-fast-artificial sequence provides simultaneous motion variations in all DoFs with large ranges of motion and with smoothly varying and periodically increasing and decreasing angular velocities, as well as direction changes. Comparing Figure 4 with Figure 3, the latter showing the Euler angle sequences and ranges of motion of real-slow and real-fast, sim-fast-artificial, sampling a range of motion of $[\pm {139}^{\circ }]$, includes the maximum range of motion reached in the captured data sequences (${108}^{\circ }$ elbow extension in real-slow), however, applies this to each rotational DoF. Moreover, comparing Table 6 with Table 2, sim-fast-artificial also reaches comparable peak angular velocity 2-norms to sim-fast. The peak accelerometer 2-norms are lower due to the artificial motion being more smooth, however, they are still considerably higher than in real-slow. Note, these accelerations depend on the assumed I2S positions. While sim-fast-artificial shares the mentioned characteristics with real-fast, as indicated above, it represents an artificial motion, which does neither resemble human motion patterns nor respect anatomical movement restrictions. However, it contains systematically sampled large ranges of motion and changing angular velocities in all DoFs, providing a challenging test case, as also visible in the supplementary video. Besides this, sim-fast-artificial also supports easy reproducibility by other researchers, without the need for sharing recorded data sequences, since the movement is based on an analytic function. Note, a similar sequence was also used in [40].

In contrast to the artificial motion in sim-fast-artificial, sim-slow and sim-fast represent re-simulations of the captured data sets real-slow and real-fast, i.e., they provide simulated IMU data for the human motions described in Section 2.6.1. As mentioned before, these movements are illustrated in the supplementary videos. For the simulation, the segment lengths, I2S calibrations and IMU trajectories were obtained from the optical system (cf. Section 2.6.1) and the IMU data was again obtained by differentiation. To ensure the IMU data to stay in feasible ranges, the IMU trajectories were low-pass filtered with a cut-off frequency of 10 Hz using a Butterworth filter of 4th order [58]. Table 2 shows the peak acceleration and angular velocity 2-norms of sim-slow and sim-fast in comparison to the respective real data sequences. Compared to sim-fast-artificial, the purpose of sim-slow and sim-fast is to provide IMU data based on human motion including respective limited joint DoFs and motion ranges, and, as already mentioned in Section 2.6.1, to provide slow and fast motion separately, in order to enable performance comparison of the different sensor fusion methods in relation to motion agility. Compared to real-slow and real-fast, the purpose of sim-slow and sim-fast is to provide IMU data for moderate and fast human motion with consistent ground truth and without errors, such as those described in Section 2.6.2. Moreover, the re-simulated data sequences enable comparison of the performances of the proposed sensor fusion methods between simulated and real data of the same motion.

Based on the above described motion sequences the influences of I2S calibration and segment length errors on the proposed sensor fusion methods were studied by systematically introducing respective model calibration errors when simulating the IMU data of the middle segment and comparing the estimated segment orientations with the ground truth data (cf. Section 2.6.4). The middle segment was chosen in order to enable assessing the error propagation behavior into neighboring segments. To this end, five calibration error types were introduced:

• I2S position errors: $\Delta p_{I2S}^{along}\in [-0.2\mathrm{m},0.2\mathrm{m}]$, i.e., position changes along the segment axis.
• I2S position errors: $\Delta p_{I2S}^{out}\in [-0.2\mathrm{m},0.2\mathrm{m}]$, i.e., position changes perpendicular to the segment axis.
• Segment length variations: $\parallel \Delta S_1\parallel \in [-0.2\mathrm{m},0.2\mathrm{m}]$.
• I2S orientation errors: $\Delta q_{I2S}^{along}\in [-{30}^{\circ },{30}^{\circ }]$ along the bone, i.e., rotations around the segment axis associated to the IMU.
• I2S orientation errors: $\Delta q_{I2S}^{out}\in [-{30}^{\circ },{30}^{\circ }]$ out of bone, i.e., rotations around the IMU axis perpendicular to the surface of the associated segment.

The ranges for the orientation errors are based on the findings in [22] concerning the trueness (deviation from reference) of established calibration methods. The reported maximum error of ${26}^{\circ }$ was slightly increased in order to account for the fact that larger errors are to be expected when calibration movements are performed by a subject autonomously. For obtaining the I2S positions and segment lengths, as described in Section 1.2, a variety of methods providing different and not well documented levels of reliability are in use. Hence, for this study, we determined the maximum I2S position error that could appear for the test subject by assuming the IMU to sit inside the joint, i.e., a $0.2$ m shift along the forearm segment. We then used this value symmetrically for all translational error ranges, in order to be able to compare the influences of the different error types also relative to each other using a subsequently described normalized error measure (cf. Section 2.6.4). Note that these error ranges likely include the majority of errors that could potentially appear, e.g., due to varying body shapes and structures, heavy clothing or also failure of self-calibration methods, while they allow to compare the robustness of the different sensor fusion methods also under extreme conditions. Note, all considered model calibration errors were introduced separately in order to assess their influence in isolation.

The above error ranges were sampled with steps of $2^{\circ }$ and $0.02$ m, respectively. In total, this summed up to 120 evaluations for each sensor fusion method.

2.6.4. Error Measures

The performances of the different sensor fusion methods were evaluated based on the angular deviations of the estimated segment orientations, $q_{i,t}^{GS}$, from the reference segment orientations, ${\overline{q}}_{i,t}^{GI}$. For the real data scenario, the latter were calculated from the optical data according to Equation (16). For the simulation scenario, the segment orientations were available from the given motion sequence and can be considered as error-free ground truth. The angular deviations for each segment i at each time step t were calculated as [59]:

$$E_{i,t}^{rel}=2\mathrm{acos}{({\overline{q}}_{i,t}^{GS}\odot q_{i,t}^{SG})}_w$$

(21)

Here, ${(·)}_w$ denotes the scalar part of the quaternion. The angular error statistics in terms of mean, standard deviation and maximum orientation errors over K segments and N time steps of a data sequence were calculated as:

$$\begin{array}{ccc}\hfill E^{mean}& =& \frac{1}{K}\sum _{i=0}^{K-1}\frac{1}{N}\sum _{t=0}^{N-1}{E}_{i,t}^{rel}\hfill \end{array}$$

(22a)

$$\begin{array}{ccc}\hfill E^{sd}& =& \frac{1}{K}\sum _{i=0}^{K-1}\sqrt{\frac{1}{N-1}{\sum }_{t=0}^{N-1}{(E^{mean}-E_{i,t}^{rel})}^2}\hfill \end{array}$$

(22b)

$$\begin{array}{ccc}\hfill E^{max}& =& \underset{i=0,\dots ,K-1}{max}\underset{t=0,\dots ,N-1}{max}{E}_{i,t}^{rel}\hfill \end{array}$$

(22c)

By setting $K=1$, the angular error statistics are calculated for each segment separately ($K=3$ considers all segments used in the evaluation). All angular errors are given in degrees.

Additionally to the above measures of performance, a normalized range error, $E_{\text{tracker}}^{norm}$, was also defined based on the mean angular errors, for quantifying the influence of model calibration errors on the different sensor fusion methods and simplifying the performance comparison. The normalized range error was defined as:

$$\begin{array}{c}\hfill E_{\text{tracker}}^{norm}=\frac{{\sum }_{i=0}^{N_{var}-1}{E}_{\text{tracker}}^{mean}}{{max}_{j=0,\dots ,N_{test}-1}{\sum }_{i=0}^{N_{var}-1}{E}_{ref}^{mean}}\end{array}$$

(23)

Here, $N_{var}$ is the number of samples in the considered calibration error range, i.e., $N_{var}=30$ for I2S orientation errors and $N_{var}=20$ for I2S position and segment length errors. Moreover, $N_{test}$ denotes the number of tests considered ($N_{test}=6$ in the evaluation, i.e., three sequences, sim-fast-artificial, sim-slow, sim-fast, and two test conditions, w/mag, w/o mag). In Equation (23), the denominator scales the sum of mean errors, such that $E_{\text{tracker}}^{norm}$ is formulated relative to a reference $E_{\text{ref}}^{mean}$. For the latter, the Chaintracker was chosen, since this was considered as baseline for our work. Hence, for $E_{\text{tracker}}^{norm}=1$, the performance of the respective method is similar to the Chaintracker. Lower values indicate a better performance, higher errors indicate a worse performance.

3. Results

3.1. Tracking Performances on Real Data

This section presents the performance results of the different sensor fusion methods on real data, taking into account the inaccuracies as described in Section 2.6.1. The per segment angular error statistics (cf. Equation (22)) for the four sensor fusion methods (Chaintracker, Quattracker segment, Quattracker IMU, Optitracker) on real-slow and real-fast under both test conditions (w/mag, w/o mag) are summarized in

Table 7. Real data scenario: mean (std; max) angular errors for each segment. Note, the color represents a linear interpolation of the mean error over all segments between red (maximum error) and green (minimum error). This helps comparing the performances of the different sensor fusion methods. Also note, w/mag refers to using the real magnetometer measurements, w/sim. mag refers to using the simulated magnetometer measurements (cf. Section 2.6.2) and w/o mag refers to dropping the magnetometer information.

. The following general tendencies could be observed: The Optitracker performed best w.r.t. mean errors, standard deviations and maximum errors on both data sequences and test conditions. The maximum error was $6.{78}^{\circ }$ and appeared at the forearm segment on real-fast. The mean errors were below $2.4^{\circ }$, and the standard deviations were below $1.4^{\circ }$ in all cases. Note the comparably low weight of the magnetometer measurements as used in the Optitracker (cf. Appendix D).

Compared to the Optitracker, the error margins of all EKF-based methods were considerably higher, if the real magnetometer measurements were used (Table 7, data rows 1–3 and 10–12, cf. Table 4). However, the error margins, in particular of the Quattracker IMU, were only slightly higher or comparable on real-slow, if the simulated magnetic field measurements (sim. mag) were used (data row 4–6). These error margins increased for all EKF-based methods on real-fast, although the mean errors remained below $4.3^{\circ }$. However, the maximum errors differed considerably. The Quattracker segment showed a maximum error of $14.{94}^{\circ }$ with the simulated magnetic fields and the Chaintracker showed a maximum error of $23.{52}^{\circ }$ with the real magnetometer measurements on real-fast. Moreover, over all segments, the Chaintracker had consistently (on both data sequences, w/mag and w/sim. mag) the highest mean errors and standard deviations, followed by the errors of the Quattracker segment, while the Quattracker IMU performed best among the EKF-based methods, indicating a higher robustness and accuracy for the IMU centered approaches in contrast to the segment centered ones. Also over all segments, all EKF-based methods showed a mean and maximum error increase on both data sequences, if the (undisturbed) magnetometer information was dropped. On sim-slow, mean errors, standard deviations and maximum errors increased on per segment basis. Note, in particular the maximum errors of the EKF-based methods increased on real-slow without using magnetometer information, but comparably less on real-fast. The Optitracker, in contrast, was not affected by dropping the magnetometer information, neither on real-slow nor on real-fast. In summary, the following indications could be observed: Under beneficial conditions (moderate motion, no magnetic disturbances), the EKF-based methods, in particular the Quattracker IMU, performed comparable or only slightly worse than the Optitracker. However, under fast motion changes (as present in real-fast) or magnetic disturbances, or when no magnetometer information was used, the performances of the EKF-based methods deteriorated, while the Optitracker was not considerably affected. Again, note the comparably low weight of the magnetometer measurements (cf. Appendix D). The higher robustness and accuracy of the Optitracker, however, comes at the cost of a considerable increase in computational complexity (cf. Table 1).

3.2. Tracking Performances on Simulated Data with Model Calibration Errors

This section presents the performance results of the different sensor fusion methods on the simulated data sequences with systematically introduced model calibration errors. In order to evaluate the influence of these errors on the segment orientation estimation accuracy in isolation, perfect IMU data was here considered (cf. Section 3.3 for an evaluation with noisy data). The normalized range errors for the four sensor fusion methods on the three sequences (sim-fast-artificial, sim-slow, sim-fast) under both test conditions (w/mag, w/o mag) and for the five calibration error types ($\Delta p_{I2S}^{along}$,$\Delta p_{I2S}^{out}$,$\parallel \Delta S_1\parallel$,$\Delta q_{I2S}^{along}$,$\Delta q_{I2S}^{out}$) are summarized in

Table 8. Simulation scenario: the normalized range error $E_{quot}^{tracker}$ (cf. Equation (23)) is shown for the different simulated model calibration errors and sensor fusion methods. The latter have the following shortcuts: Chaintracker (chain), Quattracker segment (Seg. q.), Quattracker IMU (IMU q.), Optitracker (opti.). Note, for each table separately, the color represents a linear interpolation of the error from red (maximum error) to green (minimum error).

. Moreover, the per segment mean angular errors for sim-fast are illustrated in Figure 5, Figure 6 and Figure 7. Here, sim-fast was chosen for a more detailed presentation, since it was the most challenging human motion considered. The results were analyzed w.r.t. the following aspects: (1) severity of angular error increases and error propagation from the affected middle segment $S_1$ to the previous ($S_0$) and subsequent ($S_2$) segments in relation to the error type; (2) aspect (1) w/o using magnetometers; (3) aspect (1) in relation to the motion agility, i.e., slow and fast motion and motion changes as available in sim-slow, sim-fast and sim-fast-artificial.

First, the results with magnetometers are described (Table 8a, solid plots in Figure 5, Figure 6 and Figure 7).

In Table 8a, the highest normalized range error of $0.51$ was observed for the Chaintracker on sim-fast-artificial with error type $\Delta p_{I2S}^{out}$. Note that this was due to extremely high mean errors when applying an out-of-bone shift of above $0.18$ m. This might exceed expected out-of-bone position errors, which, however, could be better handled by the other methods.

Besides this effect, the table shows, that I2S orientation errors had a higher influence on the estimation accuracy of the segment orientations than I2S position and segment length errors. This held for all data sequences and methods (lower errors in data rows 1–3 compared to 4–5). Moreover, all EKF-based methods showed an error increase from sim-slow to sim-fast and sim-fast-artificial, for all error types (data columns 5–8 compared to 9–12 and 1–4), while the Optitracker was not affected by the motion agility. In the majority of tests (combinations of data sequences and error types), the Chaintracker showed the highest errors, followed by the Quattracker segment, followed by the Quattracker IMU.

In addition to Table 8, Figure 5 illustrates the per segment distribution of mean angular errors for both orientation error types, $\Delta q_{I2S}^{along}$ and $\Delta q_{I2S}^{out}$, on sim-fast. As can be seen in Figure 5b,e, I2S orientation errors propagated at least linearly into the affected segment $S_1$. This held for both error types, $\Delta q_{I2S}^{along}$ and $\Delta q_{I2S}^{out}$, for all methods, and also for all data sequences (i.e., also for sim-slow and sim-fast-artificial, which are not shown in detail here). Note, a direct (linear) propagation is here expected, since all joints are modeled with three DoFs, without any explicit correction strategies, such as joint constraints. However, the different methods varied in the amount of error that was propagated into the neighboring segments (Figure 5a,c,d,f). While the Optitracker showed no noticeable error propagation (maximum mean error in $S_0$, $S_2$ on sim-fast: $0.{03}^{\circ }$), the EKF-based methods showed different behaviors. Regarding $\Delta q_{I2S}^{along}$, $S_0$ and $S_2$ showed errors that were rather constant over the whole range at different levels for all EKF-based methods. Regarding $\Delta q_{I2S}^{out}$, a real error propagation could be observed for all EKF-based methods, i.e., the mean errors in $S_0$ and $S_1$ varied depending on the introduced calibration errors. Comparable propagation behaviors were also observed for sim-slow and sim-fast-artificial, but with different error levels.

As already indicated above, for all EKF-based methods, the influence of I2S position errors on the segment estimation accuracy increased with higher motion agility (Table 8, data rows 1–3). In addition to the normalized range errors in the table, Figure 6 illustrates the per segment distribution of mean angular errors for both error types, $\Delta p_{I2S}^{along}$ and $\Delta p_{I2S}^{out}$, on sim-fast. Here, for all EKF-based methods, a clear influence on the affected segment (Figure 6b,e), as well as, a clear propagation into the neighboring segments (Figure 6a,c,d,f) could be observed. Note, while this observation held for both fast sequences, sim-fast and also sim-fast-artificial, the type of I2S position error, which caused the most severe normalized range error, was sequence dependent. While on sim-fast, $\Delta p_{I2S}^{along}$ caused the highest errors for all EKF-based methods, on sim-fast-artificial, $\Delta p_{I2S}^{out}$ caused the highest errors (Table 8, data rows 1–3, data columns 1–3 and 9–11). In contrast to the EKF-based methods, the Optitracker showed no noticeable error increases in any of the segments (maximum mean error on sim-fast: $0.{04}^{\circ }$), neither on sim-fast, nor on sim-slow or sim-fast-artificial. It could therefore be considered invariant to I2S position errors in the tested settings.

The same tendencies as for the I2S position errors could also be observed for all EKF-based methods for segment length errors, i.e., a normalized range error increase from sim-slow to sim-fast and sim-fast-artificial (Table 8, data row 3), an error propagation on the fast sequences (shown in Figure 7a–c for sim-fast), and also a sequence dependent level of influence of this error type. On sim-fast, $\left|\right|\Delta S_1\left|\right|$ had a lower or comparable influence than/to both I2S position error types on the normalized range errors of the EKF-based methods (Table 8, data rows 1–3, data columns 9–11), and also on the per segment mean angular errors (maximum mean angular error over all segments and methods on sim-fast: $5.{30}^{\circ }$ by the Quattracker segment). In contrast, on sim-fast-artificial, $\left|\right|\Delta S_1\left|\right|$ had a higher or comparable influence than/to $\Delta p_{I2S}^{along}$ (Table 8, data rows 1–3, data columns 1–3). Again, the Optitracker showed no noticeable error increases in any of the segments (maximum mean angular error on sim-fast: $0.{04}^{\circ }$), neither on sim-fast, nor on sim-slow or sim-fast-artificial.

In summary, in this study, I2S orientation errors (in particular $\Delta q_{I2S}^{out}$) had the most severe influence on the accuracy of the estimated segment orientations, for all methods and data sequences. For the EKF-based methods, the influence of translational error types (I2S position and segment length errors) was lower, but increasing with motion agility (i.e., higher for sim-fast and sim-fast-artificial compared to sim-slow). The specific influences of the different translational error types compared to each other were consistent between EKF-based methods, however, sequence dependent (different results for sim-fast-artificial compared to sim-fast), i.e., depending on the individual motion. The Optitracker was not considerably affected by those translational model calibration errors, independently of the data sequence. Over all data sequences, the Quattracker IMU performed best among the EKF-based methods and the Optitracker performed best among all methods.

In the following, the results without using magnetometer information are described (Table 8, dashed plots in Figure 5, Figure 6 and Figure 7). As for the real scenario, the results for the Optitracker were nearly identical under both test conditions and for all data sequences. Therefore, only the results for the EKF-based methods are further discussed in more detail.

Table 8 shows an increase of the normalized range errors for all EKF-based methods, when the (undisturbed) magnetometer information was dropped. This held for all data sequences and error types and is an expected behavior. An interesting fact is, however, that the error increase due to introduced I2S position and segment length errors was equally or more severe than the error increase due to introduced I2S orientation errors, in particular for the Quattracker segment and Quattracker IMU and for all data sequences (data columns 2–3, 6–7, 10–11). Moreover, Table 8, similarly to the test condition with magnetometers, showed an error increase from sim-slow to sim-fast and sim-fast-artificial for Quattracker segment and Quattracker IMU (data columns 6–7 compared to 2–3 and 10–11), while, again, the type of translational error that caused the highest normalized range error increase was sequence dependent. On sim-fast, segment length errors resulted in even higher normalized range errors than those caused by I2S orientation errors (data columns 10–11), while on sim-fast-artificial, $\Delta p_{I2S}^{out}$ had more influence among the translational error types. The Chaintracker, in contrast to Quattracker segment and Quattracker IMU, seemed to benefit from the agile motions in sim-fast and sim-fast-artificial, for both I2S orientation error types, but in a complementary manner for the translational error types (data column 5 compared to 1 and 9).

In addition to the table, Figure 6 and Figure 7 show exemplary for sim-fast that dropping the magnetometer information not only resulted in increased mean angular errors for the affected segment, but also in more severe error propagation into neighboring segments. This held for both I2S position and segment length errors. Concerning the I2S orientation errors (Figure 5), dropping the magnetometer information resulted in comparable performance for the affected segment (Figure 5b,e), while a higher propagation could be observed to the neighboring segments (Figure 5a,c,d,f).

A further peculiarity visible in the figures is the asymmetry of most of the error distributions over the error ranges when dropping the magnetometer information. Similar effects were observed for sim-fast-artificial and sim-slow. As exemplary shown in Figure 7d–f, the errors and, in particular the asymmetries, appeared mainly in the yaw direction. This is expected, since the yaw direction is the variable, where the magnetometer measurements provide reference information. Without this correcting information, the error in the yaw direction is supposed to be highly dependent on the individual configuration and motion leading to unpredictable results.

In summary, even on short motion sequences and noise-free IMU data, dropping the magnetometer information resulted in all EKF-based methods being more severely affected by model calibration errors. Moreover, translational error types (I2S position and segment length errors) gained more influence compared to I2S orientation errors, in particular for Quattracker segment and Quattracker IMU and increasing with motion agility. As for the test condition with magnetometers, the type of translational error that caused the most severe normalized range error was sequence dependent, but consistent between EKF-based methods. The Optitracker, in contrast to the EKF-based methods, was not affected by dropping the magnetometer information.

3.3. Tracking Performances on Simulated Data without Calibration Errors

When focusing on the results with no (or small) calibration errors in Figure 5, Figure 6 and Figure 7 (i.e., x around 0 in all graphs), the previously (real data scenario, simulation scenario w/mag) observed performance ranking (Chaintracker, Quattracker segment, Quattracker IMU, Optitracker) could be confirmed also for the simulation scenario w/o magnetometers, as exemplary shown in numbers for sim-fast in

Table 9. Simulation scenario without calibration errors (sim-fast): mean (std; max) angular errors over all segments on the two test configurations (w/mag, w/o mag), both on noise-free (perfect) and noisy IMU data.

Method	Chaintracker	Quattracker segment	Quattracker IMU	Optitracker
Noise-free w/mag	1.42${}^{\circ }$ (1.40${}^{\circ }$; 8.04${}^{\circ }$)	1.19${}^{\circ }$ (1.23${}^{\circ }$; 6.83${}^{\circ }$)	0.66${}^{\circ }$ (0.64${}^{\circ }$; 3.30${}^{\circ }$)	0.01${}^{\circ }$ (0.01${}^{\circ }$; 0.06${}^{\circ }$)
Noise-free w/o mag	3.50${}^{\circ }$ (2.57${}^{\circ }$; 9.45${}^{\circ }$)	1.57${}^{\circ }$ (1.45${}^{\circ }$; 7.52${}^{\circ }$)	0.97${}^{\circ }$ (0.65${}^{\circ }$; 3.36${}^{\circ }$)	0.01${}^{\circ }$ (0.01${}^{\circ }$; 0.06${}^{\circ }$)
Noise w/mag	1.46${}^{\circ }$ (1.39${}^{\circ }$; 8.09${}^{\circ }$)	1.22${}^{\circ }$ (1.21${}^{\circ }$; 6.83${}^{\circ }$)	0.69${}^{\circ }$ (0.62${}^{\circ }$; 3.28${}^{\circ }$)	0.40${}^{\circ }$ (0.05${}^{\circ }$; 0.49${}^{\circ }$)
Noise, w/o mag	3.73${}^{\circ }$ (2.68${}^{\circ }$; 9.76${}^{\circ }$)	1.55${}^{\circ }$ (1.40${}^{\circ }$; 7.42${}^{\circ }$)	0.95${}^{\circ }$ (0.65${}^{\circ }$; 3.31${}^{\circ }$)	0.40${}^{\circ }$ (0.05${}^{\circ }$; 0.49${}^{\circ }$)

. The table provides the angular error statistics over all segments for the four sensor fusion methods on sim-fast under both test conditions (w/mag, w/o mag) and on noise-free and noisy IMU data. Hence, the table complements the results in the previous section for sim-fast by detailing the case, where no calibration errors were simulated and by providing results with noisy data. Note, comparable results as described below were obtained from sim-slow and sim-fast-artificial, but with overall lower error levels.

With noise-free IMU data (Table 9, data rows 1–2), the errors mainly depend on linearization, the degree to which the motion model fits the actual motion and the noise settings (cf. Appendix D). Concerning the former, the Chaintracker is expected to have the highest linearization errors, due to the angle parametrization, and the Optitracker is expected to have the lowest linearization errors, due to the nonlinear optimization. This is confirmed in the table. Moreover, the above performance ranking held for all test configurations and for mean errors, standard deviations and maximum errors. From the EKF-based methods, the Quattracker IMU provided the best results, with a mean error below $1^{\circ }$ for all test configurations. The Optitracker performed overall best, with a mean error below $0.4^{\circ }$ and a maximum error of $0.{49}^{\circ }$ over all test configurations.

For the EKF-based methods, noise did not have a significant influence (data rows 1–2 compared to 3–4). For the Optitracker, adding noise resulted in the error levels increasing according to the noise levels, with the above maximum error.

As previously observed, dropping the magnetometer information had no influence on the Optitracker. Also the Quattracker segment and Quattracker IMU showed an only slight error increase (mean, std, max) on this rather short sim-fast sequence, both on noise-free and noisy IMU data. Only the Chaintracker showed a more considerable error increase (mainly mean and std) when dropping the magnetometer information (data column 1). Note, that this effect might be reduced through specific tuning or by adding joint constraints (e.g., [20]). Since mainly the biomechanical model representation (kinematic chain versus free segments) and the rotation parametrization (joint angles vs. quaternions) differ between the two methods, this finding indicates that the combination of free segment model and quaternions is potentially more robust in the case of body motion tracking with unconstrained joints and without using magnetometer information.

4. Discussion and Conclusions

The goal of this work was to develop/identify sensor fusion methods, which robustly handle model calibration errors (I2S, segment length) and show potential for magnetometer-free operation. Concerning the latter, the conducted tests served mainly for excluding methods that do not cope well with missing reference information in the yaw direction, already on the rather short data sequences considered. Here, the motivation for this investigation was that potentially large model calibration errors have to be expected given today’s most widely used calibration techniques based on static poses, functional movements, manual measurements and assumptions. At the same time, magnetic disturbances have to be expected in man-made environments.

To this end, two newly developed/adapted sensor fusion methods (Quattracker segment, Quattracker IMU, Optitracker) were compared to an existing method (Chaintracker), regarding segment orientation tracking accuracy on complex moderate and fast captured and re-simulated human motion, as well as on fast artificial motion, in the presence and in the absence of model calibration errors, with and without using magnetometer information.

The findings concerning the influence of model calibration errors on the different methods, which held independently of the considered data sequences, can be summarized as follows: With (undisturbed) magnetometer information, I2S orientation errors clearly had a higher influence on segment orientation estimation accuracy than I2S position and segment length errors, which is an expected behavior. Both tested orientation error types, $\Delta q_{I2S}^{along}$ and $\Delta q_{I2S}^{out}$, resulted in at least linear error propagation into the affected segment. This held for all methods and was also expected, since no explicit compensation strategies, such as joint constraints were used (cf. [20]). For the EKF-based methods, different amounts of errors were also propagated into the orientation estimates of the neighboring segments, in particular for $\Delta q_{I2S}^{out}$, similarly for I2S position ($\Delta p_{I2S}^{along},\Delta p_{I2S}^{out}$) and segment length errors ($\parallel \Delta S_1\parallel$), increasing with motion agility. When dropping the magnetometer information, the influence of translational calibration error types (I2S position and segment length errors), relative to the influence of I2S orientation errors, increased considerably. From this, it can be concluded that I2S orientation errors represent the dominant source of the considered model calibration errors in magnetic-inertial human motion tracking (assuming undisturbed magnetometer information), while the importance of accurate I2S position and segment length calibration increases with motion agility and in the absence of reference information for the yaw direction (w/o mag). Note that in these cases all translational error types should be kept small, since their influences turned out to be motion dependent.

Further considering the results on real and simulated data without introduced calibration errors (cf. Section 3.1 and Section 3.3), it can also be stated for the newly developed EKF-based methods that dropping the magnetometer information (at least for the short time duration considered in the tests), resulted in acceptable accuracy loss, as long as no severe model calibration errors were present at the same time. The Chaintracker, however, showed the highest dependency on undisturbed magnetometer information and was overall more fragile.

The Optitracker, in contrast to all EKF-based methods, showed no considerable I2S orientation error propagation into neighboring segments, independently of the type of orientation error introduced ($\Delta q_{I2S}^{along}$, $\Delta q_{I2S}^{out}$). Moreover, it was nearly invariant w.r.t. I2S position and segment length errors, independently of motion agility and magnetometer usage. The latter was also confirmed in the real data scenario and in the simulation scenario without calibration errors. Note, the Optitracker also turned out to require no intensive tuning (cf. Appendix D).

Among the EKF-based methods, in the majority of test cases, the Chaintracker showed the highest (mean) errors (over all segments), followed by the Quattracker segment, followed by the Quattracker IMU. The Optitracker outperformed all EKF-based methods in all test cases, also showing lower standard deviations and maximum errors.

Hence, the segment centered biomechanical models (used by Chaintracker and Quattracker segment), which assume constant angular acceleration in the motion model (cf. [16] and Equation (8)), were overall outperformed by the IMU centered free segments biomechanical models (used by both, Quattracker IMU and Optitracker), while in the former case, the free segments model overall performed slightly better than the kinematic chain model. Moreover, the additional redundancy in terms of spatial and temporal estimation variables, in combination with the nonlinear optimization used by the Optitracker provided a considerable gain in robustness, in particular in the presence of disturbing effects as described above. Obviously, this comes at the cost of a significantly higher computational complexity, where, however, a parallelized version, similar to the one proposed in [32], might give remedy. Concerning the analysis above, the following aspects have to be noted: (1) the performances of the EKF-based methods could likely be improved through sequence dependent tuning, which is, however, not preferred for a practical system; (2) there exist other motion models, e.g., the decaying angular acceleration model [20], which were not investigated here.

Another important finding was that under beneficial conditions (moderate motion, no severe calibration errors, undisturbed magnetometer information), the Quattracker IMU provided the most comparable results to the Optitracker, but with considerably lower computational costs. Hence, the Optitracker might be best invested, if disturbing effects are expected in the considered application scenario, while, under beneficial conditions, the Quattracker IMU can be a computationally more efficient solution, which provides comparable performance. Based on this work, long-term magnetometer-free operation and robust handling of soft tissue artifacts remain parts of our future work.