Joint axis estimation for fast and slow movements using weighted gyroscope and acceleration constraints

TL;DR

This paper introduces a new method for joint axis estimation in inertial motion tracking that combines gyroscope and accelerometer constraints, achieving accurate results in both fast and slow movements without relying on magnetic field assumptions.

Contribution

The paper presents a novel cost function that integrates both acceleration and angular rate constraints with adaptive weighting, improving joint axis estimation accuracy across different motion speeds.

Findings

Accurate joint axis estimation achieved within five seconds for all motion types.

The method outperforms previous approaches in both fast and slow movements.

No reliance on magnetic field assumptions, using only accelerometer and gyroscope data.

Abstract

Sensor-to-segment calibration is a crucial step in inertial motion tracking. When two segments are connected by a hinge joint, for example in human knee and finger joints as well as in many robotic limbs, then the joint axis vector must be identified in the intrinsic sensor coordinate systems. There exist methods that identify these coordinates by solving an optimization problem that is based on kinematic joint constraints, which involve either the measured accelerations or the measured angular rates. In the current paper we demonstrate that using only one of these constraints leads to inaccurate estimates at either fast or slow motions. We propose a novel method based on a cost function that combines both constraints. The restrictive assumption of a homogeneous magnetic field is avoided by using only accelerometer and gyroscope readings. To combine the advantages of both sensor types, the residual weights are adjusted automatically based on the estimated signal variances and a nonlinear weighting of the acceleration norm difference. The method is evaluated using real data from nine different motions of an upper limb exoskeleton. Results show that, unlike previous approaches, the proposed method yields accurate joint axis estimation after only five seconds for all fast and slow motions.

Tables

Table 1. TABLE I: MAD and SAD metrics from the methods using the non-weighted residuals, e ω subscript 𝑒 𝜔 e_{\omega} and e a subscript 𝑒 𝑎 e_{a} by themselves and combined without and with the proposed weighting variables. The lowest MAD values and the MAD values within 0.5 ​ ° 0.5° are highlighted in bold for each axis–movement speed combination.

		Free axis		Vertical axis		Horizontal axis
		${\hat{j}}_1$	${\hat{j}}_2$	${\hat{j}}_1$	${\hat{j}}_2$	${\hat{j}}_1$	${\hat{j}}_2$
Speed	Residuals	MAD(SAD) [∘]	MAD(SAD) [∘]	MAD(SAD) [∘]	MAD(SAD) [∘]	MAD(SAD) [∘]	MAD(SAD) [∘]
Fast	$e_{\omega }$	2.2(1.6)	81.9(87.4)	1.4(1.7)	90.6(87.2)	8.2(5.4)	90.3(85.4)
	$e_a$	5.9(3.7)	5.9(3.3)	9.7(6.7)	10.3(6.5)	50.6(47.3)	63.4(52.0)
	$e_{\omega }+e_a$	7.6(6.7)	7.6(6.8)	4.9(3.7)	5.2(3.7)	31.6(31.2)	70.9(74.0)
	Weighted	0.7(0.4)	0.5(0.3)	0.7(0.8)	1.5(1.3)	8.3(7.4)	90.8(85.7)
Slow	$e_{\omega }$	0.8(0.4)	75.2(87.1)	2.0(1.2)	68.6(85.8)	3.4(3.0)	90.9(88.5)
	$e_a$	6.7(13.4)	6.5(14.2)	6.3(7.6)	6.9(7.4)	51.5(46.6)	58.3(47.5)
	$e_{\omega }+e_a$	3.3(11.5)	4.1(12.2)	4.4(2.4)	3.1(1.9)	48.2(35.0)	77.9(52.1)
	Weighted	0.7(0.4)	1.4(1.0)	1.0(0.5)	1.6(1.0)	2.4(1.9)	90.9(88.5)
Mixed	$e_{\omega }$	1.4(1.0)	76.7(87.0)	2.5(1.8)	79.8(84.6)	5.2(3.1)	90.6(86.9)
	$e_a$	26.2(17.4)	24.8(16.1)	17.7(11.5)	18.0(12.1)	39.4(24.9)	73.3(53.0)
	$e_{\omega }+e_a$	10.0(6.2)	7.8(5.4)	7.6(5.3)	9.5(6.5)	27.3(22.9)	82.1(70.0)
	Weighted	0.9(0.6)	1.6(1.2)	0.9(0.8)	2.1(1.5)	4.8(3.0)	90.7(87.3)

Equations

$$a_i^{S_i}$$

$$a_0^G$$

$$a_i$$

$$K(\omega ,\dot{\omega })=\left[\begin{array}{ccc}-{\omega }_y^2-{\omega }_z^2& {\omega }_x{\omega }_y-{\dot{\omega }}_z& {\omega }_x{\omega }_z+{\dot{\omega }}_y\\ {\omega }_x{\omega }_y+{\dot{\omega }}_z& -{\omega }_x^2-{\omega }_z^2& {\omega }_y{\omega }_z-{\dot{\omega }}_x\\ {\omega }_x{\omega }_z-{\dot{\omega }}_y& {\omega }_y{\omega }_z+{\dot{\omega }}_x& -{\omega }_x^2-{\omega }_y^2\end{array}\right],$$

$$R_1{j}_1$$

$${\omega }_i$$

$${\omega }_{j_i}$$

$${\omega }_{j_i^{\perp }}$$

$$\Vert {\omega }_{j_1^{\perp }}\Vert$$

$$\Vert {\omega }_i-j_i^{\top }{\omega }_i{j}_i\Vert$$

$$\Vert {\omega }_1\times j_1\Vert -\Vert {\omega }_2\times j_2\Vert =0,$$

$$j_i^{\top }{a}_i{j}_i$$

$$\begin{array}{rl}{R}_1{j}_1{j}_1^{\top }{a}_0^{S_1}& =R_2{j}_2{j}_2^{\top }{a}_0^{S_2}⟹\\ j_1^{\top }{a}_0^{S_1}& =j_2^{\top }{a}_0^{S_2}.\end{array}$$

$$\begin{array}{rl}{j}_1^{\top }{a}_1-j_2^{\top }{a}_2& =j_1^{\top }{a}_0^{S_1}-j_2^{\top }{a}_0^{S_2}+j_1^{\top }{K}_1{r}_1-j_2^{\top }{K}_2{r}_2\\ & =j_1^{\top }{K}_1{r}_1-j_2^{\top }{K}_2{r}_2,\end{array}$$

$$j_1^{\top }{a}_1-j_2^{\top }{a}_2$$

$$x$$

$$j_i(x)$$

$$\hat{x}$$

$$V(x)$$

$$e_{\omega }(k,x)$$

$$e_a(k,x)$$

$$\frac{\partial e_{\omega }(k,x)}{\partial x}$$

$$\frac{\partial e_{\omega }(k,x)}{\partial j}$$

$$\frac{\partial (\Vert {\omega }_i(t_k)\times j_i(x)\Vert )}{\partial j_i}$$

$$\frac{\partial e_a(k,x)}{\partial x}$$

$$\frac{\partial e_a(k,x)}{\partial j}$$

$$\frac{\partial j}{\partial x}$$

$$\frac{\partial j_i}{\partial x_i}$$

$$\omega (t_k)$$

$$a(t_k)$$

$$\mathrm{Var}\{e_{\omega }(k,x)\}$$

$$\frac{\partial e_{\omega }(k,x)}{\partial \omega (t_k)}$$

$$\Vert (j_i\times {\omega }_i(t_k))\times j_i\Vert$$

Full text

I Introduction

Microelectromechanical systems (MEMS) inertial measurement units (IMUs) are popular sensors, used to track movement in many different applications due to their mobility, low cost and small size. Areas where inertial sensors are used include, but are not limited to, medical rehabilitation, sports, robotics and navigation [1]. A typical IMU measures 6 degrees of freedom (DOF), using a 3D accelerometer to measure linear acceleration (translational movement) and a 3D gyroscope to measure the rate of change in the sensor’s orientation (rotational movement). In theory, a 6-DOF IMU measures all the quantities needed to fully track rigid body movement. However, since the IMU measurements must be integrated to obtain orientation and position estimates, any systematic error in the signal will accumulate. Orientation and position estimates are therefore typically only accurate for short periods of time, something which is often referred to as integration drift.

To alleviate the effects of integration drift, many different techniques have been proposed. Stationary accelerometer measurements can be used as absolute measurements of the inclination of the sensor, i.e. the angle of the sensor with respect to (w.r.t.) the global vertical axis. Using an additional 3D magnetometer sensor, which measures the local magnetic field, can give an absolute heading measurement, similar to a compass [2]. A 6-DOF IMU combined with a 3D magnetometer is sometimes referred to as a 9-DOF IMU, as the magnetometer is typically soldered to the standard IMU chip. For position tracking, additional information can be derived from multiple measurements of ranges/distances to base stations at known locations. This includes Global positioning systems (GPS), and local ranging positioning systems, using for example ultrawideband ranging [3]. In many applications however, adding additional hardware to compensate for the integration drift is not feasible, since the system then becomes too obtrusive or expensive. In such situations it becomes advantageous to make use of a model of the mechanical properties of the tracked object.

Kinematic chains are often used to model human body segments or robotic manipulator systems. The kinematic chain consist of multiple rigid body segments connected by joints with 1 to 3 DOF. The DOF of a joint determines how many free variables (joint angles) that are required to fully describe the relative orientation of the adjacent segments. For kinematic chains, one can then exploit prior knowledge about the sensors’ position w.r.t. the joints and relative orientation w.r.t. the joint axes of 1-DOF and 2-DOF joints to simplify the tracking problem, or make it easier to combine the information from multiple IMUs attached to the various segments of the kinematic chain (see for example [4] and [5]).

Such prior knowledge can of course be obtained by accurately measuring the placement of each sensor. However, a more user-friendly, less restrictive and less time-consuming solution is to allow the sensors to be placed arbitrarily and then use the available measurements to identify the sensor placement. Such methods have been proposed and used for the identification of the joint center of a 3-DOF joint [6, 7, 8] or the joint axis of a 1-DOF [6, 9] and a 2-DOF joint [10]. What these methods have in common is that they make use of the kinematic constraints of the underlying kinematic chain model to relate the measured accelerations and angular velocities to the unknown quantities to be identified, i.e. they fit the kinematic model to the data.

The current paper deals with the problem of estimating the direction of the axis of rotation in a joint connecting two segments. This joint axis estimation problem is further illustrated in Figure 1. In addition to the usefulness of a kinematic model for tracking the motion of the mechanism, such estimation of the axis of rotation is important in positioning and adjusting exoskeletons [11] so that the axes of rotation of the exoskeleton and the human match.

We propose a novel method for joint axis estimation that overcomes limitations of previously proposed methods. The method in [6] uses a gyroscope-based constraint, which requires that both segments rotate fast enough to let gyroscope bias become negligible. Moreover, even in the best case, the gyroscope-based constraint yields joint axis coordinates with ambiguous sign pairing [12]. Recently, a method that combines gyroscope and accelerometer readings has been proposed [13], however it requires stationary accelerometer measurements. We overcome this limitation by using a more general acceleration-based constraint that holds when the rotational and tangential acceleration components with respect to the joint axis are negligible. To balance between gyroscope and accelerometer information for movements of arbitrary speed, we use a weighted nonlinear optimization approach. More precisely, each data sample is weighted according to how much information about the joint axis it holds. This weighting approach allows accelerometer information to be included without requiring measurements under stationary conditions, and we will demonstrate that it provides accurate estimates with correct sign pairing for both fast and slow motions.

The remainder of the paper is structured as follows; The kinematic constraints of a hinge joint system are derived in Section II. In Section III we present the nonlinear optimization method, which makes use of the weighted kinematic constraints to estimate the joint axes. The proposed method is evaluated using experimental data collected from an exoskeleton frame moved by a human. The experimental data, the evaluation method and metrics are explained in Section IV. Results are presented in Section V and then discussed in Section VI. Finally, conclusions are drawn and potential future work outlined in Section VII.

II Kinematic model

II-A Kinematic constraints of two segments in a a kinematic chain

Consider first the kinematic chain model where we have two rigid body segments connected by a joint. The joint can be either 1-, 2- or 3-DOF. Furthermore, we consider the case where each segment has one IMU rigidly attached to it in an arbitrary position and orientation. We define two Cartesian coordinate frames, referred to as the sensor frames, $S_{1}$ and $S_{2}$, that are fixed in the center of the accelerometer triad of each IMU. We also define a global coordinate frame, $G$, which is fixed with respect to the local environment. We let subscripts $i\in\{1,2$} denote quantities belonging to a specific sensor frame. We have that

[TABLE]

where $a_{i}$ are the accelerations of the sensor frames, $a_{0}$ is the acceleration of the joint center, $\omega_{i}$ and $\dot{\omega}_{i}$ are the angular velocities and angular accelerations of the sensor frames. The positions of the joint center with respect to each sensor frame are denoted by $r_{i}$, which we assume to be unknown. Each quantity in (1) are vectors in $\mathbb{R}^{3}$ since they explain arbitrary movements of a real-world object. The superscripts $S_{i}$ are used to indicate when a certain quantity is expressed in one of the two sensor frames. If we know the orientation of each sensor frame with respect to the global frame we can express the quantities in the global frame instead. For example

[TABLE]

where $R_{i}^{GS_{i}}$ are the $3\times 3$ rotation matrices that maps a vector expressed in $S_{i}$ into the global frame. The multiplication between a rotation matrix and a vector is equivalent to an orthonormal change of basis.

For convenience we shall for the remainder of this document drop the use of the superscripts except for where it’s needed. The relationship in (1) is linear in $a_{0}$ and $r_{i}$ and can equivalently be formulated as

[TABLE]

where

[TABLE]

where subscripts $x,y,z$ denote the elements of the three-dimensional vectors. For convenience of notation we will write $K_{i}=K(\omega_{i},\dot{\omega}_{i})$.

II-B Kinematic constraints of a hinge joint system

For a 1-DOF joint, the two segments can only rotate independently with respect to each other along the joint axis. We let $\|\cdot\|$ denote the Euclidean vector norm, then the joint axis is defined by the unit vector $j\in\mathbb{R}^{3},\|j\|=1$. We refer to such a joint as a hinge joint. From its definition, the joint axis vector, $j$, must satisfy

[TABLE]

meaning that the joint axes $j_{i}$ expressed in the two sensor frames coincide in the global frame, see Figure 2. We can decompose the angular velocities into one component that is parallel to the joint axis and one that is perpendicular to the joint axis

[TABLE]

Since the two segments can only rotate independently along the joint axis it follows that the perpendicular components must have the same magnitude regardless of reference frame

[TABLE]

The magnitude of the perpendicular component can also be computed from the cross product between the angular velocity and the joint axis

[TABLE]

Combining (9) and (10) we formulate the angular velocity constraint

[TABLE]

which must be satisfied by hinge joint systems.

If we now look at the projection of the accelerations onto the joint axis, from (3) we have that

[TABLE]

It then follows from (5) that

[TABLE]

By projecting the accelerations onto the joint axis and subtracting one from the other we get

[TABLE]

where we see that only the rotational components of the accelerations remain on the right hand side. It’s clear that if the radial and tangential acceleration components along the direction of the joint axis are small ($j_{i}^{\top}K_{i}r_{i}\approx 0,\,\forall{i}$) the right hand side will vanish

[TABLE]

which forms the acceleration constraint for the hinge joint system.

III Estimation method

III-A Nonlinear optimization

Since the joint axes $j_{i}$ are unit vectors, we parametrize them using spherical coordinates to implicitly enforce the unit vector constraint

[TABLE]

which then become the unknown parameters to estimate. The estimation problem can then be formulated as

[TABLE]

where $e_{\omega}(k,x)$ and $e_{a}(k,x)$ are residual terms, based on the angular velocity- and acceleration constraints given by (11) and (15)

[TABLE]

Important to note here is that we now use the IMU measurements of the kinematic variables $\omega_{i}$ and $a_{i}$, where $t_{k},\,k=1,\ldots,N$ denotes the sampling time of the $k$:th data sample used for estimation. The two scalars $w_{\omega}(k)$ and $w_{a}(k)$ are weights that are used to give more or less importance to specific samples. A strategy for selecting these weights is proposed in Section III-B.

The optimization problem (18) is a minimization problem for a sum of squared residual terms. An efficient method for solving such problems is the Gauss-Newton method [14], which uses the gradients of the residuals to estimate the Hessian of the cost function, which is then used to select the search direction of the optimization method.

The gradients of the residuals are computed in the following way

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

III-B Weighting the residuals

The weights $w_{\omega}(k)$ and $w_{a}(k)$ that appear in the residuals (20)–(21) should in some sense reflect the amount of information the $k$:th sample contain about the joint axis. In this section we propose a strategy for selecting these weights.

The first thing to consider is that the two residuals (11) and (15) contain measurements of two different physical quantities, angular velocity [$\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$] and acceleration [$\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2}$]. Because of this, we expect the variance of the residuals to have different magnitudes. To estimate the variance of the residuals we assume that the gyroscope and accelerometer measurements are random variables that belong to multivariate Gaussian distributions

[TABLE]

where $\mu_{\omega,a}$ are the means, which depend on the movement that each sensor is undergoing at time $t_{k}$ and the $\Sigma_{\omega,a}$ are covariance matrices, assumed to be diagonal and time-independent. The means are considered unknown since they depend on the instantaneous angular velocities and accelerations of the sensors. The covariance matrices, however, can be estimated as the sample covariance of a stationary sensor. Which is typically done as part of the IMU calibration.

Because the gyroscope residual is nonlinear, we approximate its variance using the delta-method, which is based on a first-order Taylor approximation of a nonlinear function of a random variable around its mean value

[TABLE]

Important to note here is that the partial derivatives should be evaluated at the mean $\mu_{\omega}(t_{k})$, but this is not known since it depends on the true angular velocity of the system. However, the partial derivatives in (34) are unit vectors since

[TABLE]

If we then choose to approximate the covariance matrix of the measurements as

[TABLE]

where $I$ is the unit matrix and $\sigma_{\omega,a}^{2}$ is the maximum sample variance among all sensor axes. We then get a worst-case estimate of the residual variance, which is independent of the unknown sensor axes

[TABLE]

The variance of the accelerometer residual is easier to compute as it is linear. If we use the same worst-case variance approximation as for the gyroscope residual, we get

[TABLE]

We use the ratio of the estimated variances to select a base weight, $w_{0}$ as

[TABLE]

which is equal to the ratio of the worst-case sample standard deviations of the accelerometers and gyroscopes. Now if we select $w_{\omega}(k)=w_{0}$, the variances of the residuals will have similar magnitudes when both kinematic constraints are valid.

However, as previously shown, the acceleration constraint implies $j_{1}^{\top}K_{1}r_{1}\approx 0$ and $j_{2}^{\top}K_{2}r_{2}\approx 0$. Since we do not assume any prior knowledge of $j_{i}$ nor $r_{i}$ we will use the difference in magnitude of the two accelerometer measurements as a penalty term to lower the accelerometer residual weights. We use the fact that

[TABLE]

to select the accelerometer residual weight as

[TABLE]

An interpretation of this weighting strategy is that we choose $(\|a_{1}(t_{k})\|-\|a_{2}(t_{k})\|)^{2}$ as the pseudo-variance of an error term that invalidates the acceleration constraint. Note that accelerations that have equal linear and rotational components will also yield $w_{a}(k)=1$ even though they violate the constraint. However, when (41) is not satisfied, we know that the acceleration constraint is also not satisfied. An alternative strategy could be to select the error term based on $K_{1}r_{1}-K_{2}r_{2}$, but since the rotational acceleration components depend on the unknown joint center positions $r_{i}$, we deemed the chosen $w_{a}(k)$ to be less restrictive.

IV Experimental evaluation

IV-A Data acquisition

Data were collected from two Xsens wireless IMUs [15] attached to different segments of an upper limb exoskeleton frame that were connected by a hinge joint, see Figure 3. The exoskeleton was controlled by a human, performing motions where both segments move and the joint undergoes repeated flexion and extension at varying speed. First, a completely arbitrary motion is performed, where the entire joint rotated freely in 3D space. Then the joint moves such that the joint axis remains vertical while both segments move in the same horizontal plane. Finally, that planar motion is performed again but with a horizontal joint axis. The first motion is chosen to show that the methods work for arbitrary motions, while the second is chosen to show that it is not required to move the joint around all degrees of freedom. The third motion is included because recent theoretical analysis shows that for such motions the joint axis coordinates are only identifiable up to sign pairing [12]. A total of nine data sets were collected:

• •

Freely moving joint axis: Fast, slow and mixed movement speed

• •

Vertical joint axis: Fast, slow and mixed movement speed

• •

Horizontal joint axis: Fast, slow and mixed movement speed

The fast data has angular velocity magnitudes up to $4.1\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$ (Mean $=$1.3\text{,}\mathrm{rad}\text{,}{\mathrm{s}}^{-1}, SD $=$0.9\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}) and the slow data up to $0.6\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}$ (Mean $=$0.09\text{,}\mathrm{rad}\text{,}{\mathrm{s}}^{-1}, SD $=$0.08\text{\,}\mathrm{rad}\text{\,}{\mathrm{s}}^{-1}). The mixed data alternates between fast and slow movement in intervals of $5\text{\,}\mathrm{s}10\text{\,}\mathrm{s}$.

Before each data set we collected $10\text{\,}\mathrm{s}$ of measurements where the sensors were stationary. This stationary data was used to calibrate the IMUs. A bias parameter was estimated as the sample mean of the stationary gyroscope measurements. This constant bias was then subtracted from each measurement of the respective following data sets. For the accelerometer a scalar gain parameter was estimated using the least squares method in order to make the average magnitude of the stationary accelerometer measurements match that of the gravitational acceleration, $g=$9.81\text{,}\mathrm{m}\text{,}{\mathrm{s}}^{-2}$$. Each accelerometer measurement was then multiplied by this gain parameter.

IV-B Evaluation method and metrics

To evaluate the use of the proposed residual weighting method with the non-weighted alternatives we perform the following steps for each specific data set:

Divide the data set into $M$ segments spread out evenly over the length of the entire data set, where each segment contains $N$ time-consecutive measurements, i.e. $t_{1}<t_{2}<\ldots<t_{N}$. 2. 2.

For each data segment, one random initial estimate $\hat{x}_{0}$ is sampled such that the initial estimates of $j_{i}$ are distributed uniformly on the unit sphere in $\mathbb{R}^{3}$. 3. 3.

Compute the joint axis estimates $\hat{j}_{i}$ using all of the different methods to be evaluated. 4. 4.

After all $M$ estimation runs have finished, compute evaluation metrics.

The metrics used for evaluation are the mean angular deviation (MAD) and standard angular deviations (SAD) [13]. The angular deviation (AD) is defined as the angle between two estimated joint axis vectors

[TABLE]

where subscripts $m,n\in\{1,\ldots,M\}$ are used to denote specific estimates in the set of all estimated joint axes. The AD is then used to compute the MAD and SAD for $\hat{j}_{i}$ for each method to be evaluated. The MAD is computed as the sample mean of all different ADs in the set of estimated joint axes

[TABLE]

and the SAD is computed as the sample standard deviation of the same set of ADs. Since we use a rigid mechanical exoskeleton frame with orthogonal edges for evaluation it is possible to place the sensors such that we have fairly accurate knowledge about the true joint axes $j_{i}$. The sensors were placed such that one of the edges of the sensor casing was as closely aligned with the joint axis as possible. Since the edges of the sensor casing are approximately aligned with the true sensor axes, the true joint axis should have an element close to $\pm 1$ at the element corresponding to the aligned sensor axis and the other two elements should be close to zero. Using this we compute the angular error as $\text{AD}(\hat{j}_{i,m},j_{i})$ and compute the means and standard deviations as a metric for how consistent the estimates are with respect to the true joint axis.

Note that any joint axis can be described equally by a vector $j$ and its negative counterpart $-j$. Hence we consider both $(j_{1},j_{2})$ and $(-j_{1},-j_{2})$ as valid joint axis coordinates (correct sign pairing). However, it is important to note that $(-j_{1},j_{2})$, $(j_{1},-j_{2})$ are non-valid estimates (wrong sign pairing), as also discussed in [12]. For example, joint axis coordinates with wrong sign pairing lead to completely wrong results when used for calculating joint angles [5].

The AD metric should reflect the correct identification of the axis direction as well as the correct sign pairing. To assure this, we first compute $\text{AD}(\pm\hat{j}_{1,m},j_{1})$ (with respect to the true joint axis) and select the sign which yields the lowest AD. If we then select a negative sign for $\hat{j}_{1,m}$ we do the same for $\hat{j}_{2,m}$ to maintain the original sign pairing. Doing this will cause the estimates of $j_{1}$ to be as consistent as possible, whereas the estimates for $j_{2}$ will only be consistent if the sign pairing is correct.

V Results

The evaluations were performed for all $9$ different data sets, where each data set was divided into $M=100$ data segments with $N=500$ consecutive measurements in each segment. For the IMUs we used, this corresponded to $5\text{\,}\mathrm{s}$ of consecutive measured motion. The methods being compared use the gyroscope and accelerometer residuals alone, combined without weighting and combined with the weighting proposed in Section III-B. The resulting MAD and SAD values are summarized for all methods and data sets in Table I. The angular error metrics, which compare the joint axis estimates to the joint axis to which the sensors were aligned, are summarized in Figure 4.

VI Discussion

The compared methods were evaluated by their achieved MAD/SAD and angular errors. Lower MAD/SAD values indicate that the estimates are more consistent across the entire data set whereas a lower angular error indicates that the estimates are consistent with the joint axis to which the sensors were aligned. Even though the edges of the sensor casing was carefully aligned with the edges of the exoskeleton frame and the joint axis, we still expect some uncertainty in the alignment. Partially due to the placement itself and partially due to the fact that the true sensor axes may not be perfectly aligned with the sensor casing. A conservative estimate of this uncertainty is expected to be around [math]. Therefore, we cannot draw conclusions based solely on angular errors if their magnitudes are close to that of the uncertainty. Angular errors larger than [math] were deemed too inconsistent with the sensor placement, which is why Figure 4 only shows errors less than that.

VI-A Performance evaluation of the compared methods

We first compare the results for the free and vertical axis data. The proposed method with the weighted gyroscope and accelerometer residuals had the best performance out of all compared methods. This is indicated by the bold faced MAD values in Table I and also by the lowest angular errors in Figure 4. This shows that the proposed method is both the most consistent method which also achieves correct sign pairing for the free and vertical axis data.

The gyroscope-only method performed similarly to the proposed method for $\hat{j}_{1}$, which were chosen as the minimizer of $AD(\pm\hat{j}_{1},j_{1})$. However, because the sign pairing is ambiguous for the gyroscope-only method [12] we see very large MAD/SAD and angular errors for $\hat{j}_{2}$. The metrics for $\hat{j}_{1}$ demonstrate the ideal performance of the method when we have prior knowledge about the sign pairing of the joint axes, for example by restricting the sensor attachment to a known surface, and the sensor facing a predetermined side of the segment.

The accelerometer-only and the combined non-weighted methods performed significantly worse in terms of MAD/SAD and angular errors, but achieves the correct sign pairing as indicated by the lower angular errors for $\hat{j}_{2}$ for the slow movement. Interesting to note that in some cases, these methods yielded relatively low MAD/SAD values, for examples the Free axis–Fast data in Table I. However, the angular errors for the same data were significantly larger as shown in Figure 4. This indicates that these methods have larger bias for fast and mixed movement speeds, which can be explained by the fact that the approximation used to define the acceleration constraint in (14)–(15) is not valid for fast rotations along the joint axis. Combining the accelerometer and gyroscope yielded better results than accelerometer-only method, but the non-weighted combination was significantly worse than the weighted combination.

For the horizontal axis data we see that the proposed method and the gyroscope-only method performed similarly and that the accelerometer-only and combined non-weighted methods performed significantly worse. As predicted by the results concerning identifiability of the joint axis in [12], the accelerometer only allows identifiability up to sign pairing in the case where the joint axis is horizontal. A possible explanation for the worse results for $\hat{j}_{1}$ could be that the vertical acceleration component is much larger than the horizontal acceleration component (gravitational acceleration $g$ compared to the horizontal acceleration induced by flexion/extension of the joint). This might cause the scalar product between the acceleration vectors and horizontal vectors other than the joint axis, to be close to zero, causing ambiguity in the determination of a global minimum of the cost function.

The results have shown that, out of the compared methods, the proposed method is the only one that consistently achieves the lowest angular error for both $j_{1}$ and $j_{2}$ for all motions considered. Even in the horizontal axis case, where the accelerometer do not provide information about the sign pairing, the proposed method has at least similar performance to that of the gyroscope-only method. However, these results should encourage the user to not restrict the joint axis to the horizontal plane while collecting the estimation data, when it is possible.

VI-B Choice of weighting variables

The proposed weighting serves two purposes. First, the gyroscope weight $w_{\omega}(k)$ scales the gyroscope residuals to make the variance match that of the accelerometer residual. For our sensors we identified the standard deviations of the measurement noise to be $\sigma_{\omega}=$0.0050\text{,}\mathrm{rad}\text{,}{\mathrm{s}}^{-1} for the gyroscope and $\sigma_{a}=$0.0346\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-2} for the accelerometer. Which means that $w_{\omega}(k)=w_{0}=\frac{\sigma_{a}}{\sigma_{\omega}}=$6.92\text{,}\mathrm{m}\text{,}{\mathrm{rad}}^{-1}\text{,}{\mathrm{s}}^{-1}$$. Therefore, the gyroscope residuals are given more weight after this constant scaling, which prevents the accelerometer residuals from having a significantly larger impact on the value of the cost function. The second purpose is to give less weight to accelerometer residuals that violate the approximation of the constraint $\eqref{eq:acc_constrant_nonzero}$. Therefore, we chose to include a penalty term that depend on the difference in magnitude of the two accelerometer measurements in (42). However, a potential weakness of the proposed weighting strategy is that some measurements that violate the acceleration constraint might still get a weight $w_{a}(k)=1$, if the rotational acceleration components of both sensors are the same. Another approach that was considered was to select a penalty term to match the order of magnitude of $j_{1}^{\top}K_{1}r_{1}-j_{2}^{\top}K_{2}r_{2}$, since this is the actual approximation error of the acceleration constraint. However, this approach assumes prior knowledge of the position of each sensor with respect to the joint axis (the vectors $r_{i}$). It is known that identification of $r_{i}$ is ambiguous for hinge-joints [8], since in that case any vector that points from the origin of the sensor frame to any point on the infinitely extended joint axis satisfies the kinematic constraints in (1). Therefore, the chosen approach was deemed less restrictive. The proposed weighting is free from tuning parameters that need to be manually selected by the user. It depends only on the noise standard deviations, which are typically computed as part of IMU calibration, and the measurements themselves.

VII Conclusions and future work

A novel method for estimating the joint axis of a hinge joint in a kinematic chain has been proposed. The method uses a weighted combination of accelerometer and gyroscope readings obtained from two 6-DOF IMUs, one attached to each segment of the hinge joint. The proposed method was compared to methods which use gyroscope- or accelerometer information alone, or combined but non-weighed. The methods were evaluated using real data collected from sensors strapped to a rigid upper limb exoskeleton frame. The data consisted of nine different data sets, with different combinations of fast, slow and mixed movement speed under conditions where the joint axis was either free to move or constrained to a horizontal or vertical orientation. In contrast to all other methods, the proposed method yields accurate joint axis estimates for all movement speeds, under the free- and vertical joint axis conditions. For the horizontal joint axis condition, the proposed method performs similarly to the gyroscope-only method when sign pairing cannot be determined by the accelerometer.

The proposed method can be used in any system with approximate hinge joints without the need for accurate sensor placement or precise calibration motions. The inclusion of the weighted accelerometer information overcomes many limitations of the previously proposed methods. Therefore, the proposed method is an important step towards full plug-and-play motion tracking, where arbitrarily placed sensors calibrate themselves with minimal interference required by the user.

Not included in the present study is an evaluation on human data, where soft tissue artifacts may be present and cause the rigid body assumption to be violated. We have recently collected data from human movement of knee- and finger joints, which will be used in further evaluation of the method. Another interesting extension of the method would be to combine joint axis estimation with joint center estimation.