1 Introduction

Traditionally, biomechanical studies employ optical motion tracking systems for the determination of position and orientation in a local room-based coordinate system. This constrains the experiments to the calibrated volume of the camera system, although the cameras may move. As an alternative to the optical motion tracking system, several research groups propose the use of inertial sensors (accelerometers and gyroscopes) as an alternative. Besides the many advantages of these sensors compared to optical measurement systems, the inherent drift due to the unavoidable integration over time of sensor signals to obtain position and orientation introduces large errors. Moreover, it is not possible to estimate positions of inertial sensor modules with respect to each other. Although the integration drift can be reduced by using suitable estimation algorithms [7, 10, 13, 23], by using the orientations of individual body segments in a linked segment model [11, 24], or by applying suitable initial and final conditions and a limited integration time [17, 20], a stable and robust solution to estimate relative positions of sensors with respect to each other is required.

The estimation of relative positions between body segments, preferably in an ambulatory environment, is important in many applications. An example is the relation between the center of mass and the center of pressure for balance assessment. In a previous study, we proposed a method to estimate this relation by using shoes instrumented with force/moment sensors and inertial sensors [19]. The results were promising, but the relative position between the feet as well as the vertical distance between center of mass and center of pressure could not be assessed. Other examples are the estimation of relative positions in virtual reality applications [5], or to quantify mechanical loading [6].

A solution for the estimation of relative body positions is to fuse an inertial sensor system with a magnetic tracking system. An advantage of fusion with magnetic tracking [5] compared to other tracking solutions such as optical [12], acoustical [8], ultra-wide band [22], or GPS [21] is that it does not suffer from loss or degradation of the signal which can occur especially in indoor environments, or line-of-sight problems since the human body is transparent for the magnetic field applied. Roetenberg et al. [15] proposed an ambulatory magnetic position and orientation measurement system, which was fused with an inertial sensor system using a complementary Kalman filter structure in a successive study [14]. Although the system is able to estimate relative body positions and orientations accurately using a measurement system worn on the body, some important aspects of the system need to be improved. First, the estimation algorithm is based on the dipole approximation of the source and requires all three coils to be mounted orthogonally and share the same origin. Second, all source coils need to be actuated every update with a fixed update rate. Third, the fusion filter first estimates the position and orientation from the magnetic measurements which are then fed into the fusion filter as a measurement input. This means, the stochastic characteristics of the magnetic measurement system are not propagated through the fusion filter.

In order to assess these aspects, this study proposes and evaluates an alternative algorithm for relative position and orientation estimation that does not depend on a fixed coil configuration or dipole approximation, allows an optimal choice of actuation parameters, and uses the actual magnetic field that is measured as an input to the fusion filter. A complementary Kalman filter structure is presented that uses a measurement model which has been presented in a previous study [18]. The filter predicts the position and orientation based on the signals measured by the accelerometer and gyroscope of the inertial sensor. If the uncertainty associated with the relative position or orientation exceeds a predefined threshold, the system decides to perform a magnetic actuation. Moreover, only the coil that delivers most information is actuated. This way the system achieves high accuracy at relatively low energy consumption.

2 Methods

2.1 Relative position and orientation determination

Change of position and orientation can be estimated by integration of acceleration and angular velocity signals obtained from inertial sensor modules. In this study, inertial sensing is fused with a magnetic measurement system to estimate the relative position and orientation of the sensor with respect to the magnetic source, \(\user2{p}^{c}_{s}\) and \(\user2{R}^{c}_{s},\) respectively. Figure 1 shows the measurement system used to estimate the relative position and orientation of an Inertial and Magnetic Measurement System (IMMS) with respect to the magnetic source. The IMMS contains a 3D accelerometer, a 3D gyroscope, and a 3D magnetometer. A schematic overview of a configuration with a coil around the z-axis is shown in Fig. 2. The global frame is denoted by \(\Uppsi_{g},\) the magnetic source frame by \(\Uppsi_{c},\) the coil frame by \(\Uppsi_{cz},\) the sensor frame by \(\Uppsi_{s},\) and the estimated sensor frame by \(\hat{\Uppsi}_{s}.\) It should be noted that the coil frame \(\Uppsi_{cz}\) is rigidly connected to magnetic source frame \(\Uppsi_{c}.\) The reason to include both is that the magnetic source can have multiple coils attached to it with individual relative positions and orientations with respect to the magnetic source. Both the movement of the source and sensor should be estimated, since the source and sensor can move independently. Orientation is estimated by integration of angular velocity using the following differential equation [1]:

$$ \dot{\user2{R}}^{g}_{s} = {\user2{R}}^{g}_{s}\tilde{\varvec{\omega}}^{s,g}_{s}, $$
(1)

where \(\user2{R}^{g}_{s}\) denotes the rotation matrix describing the orientation of the sensor frame \(\Uppsi_{s}\) with respect to the global frame \(\Uppsi_{g}.\) The columns of \(\user2{R}^{g}_{s}\) are the unit axes of frame \(\Uppsi_{s}\) expressed in frame \(\Uppsi_{g}: {\user2{R}}^{g}_{s} =\left({\user2{X}}^{g}_{s} \;{\user2{Y}}^{g}_{s}\; {\user2{Z}}^{g}_{s}\right).\) The angular velocity \({\varvec{\omega}}^{s,g}_{s}\) of frame \(\Uppsi_{s}\) with respect to \(\Uppsi_{g},\) expressed in \(\Uppsi_{s}\) as measured by the inertial sensor is represented in skew-symmetric matrix form, indicated by the tilde operator (∼):

$$ \tilde{\varvec{\omega}}^{s,g}_{s} =\left( \begin{array}{lll} 0 & - \omega _{z} & \omega _{y} \\ \omega _{z} & 0 & - \omega _{x} \\ - \omega _{y} & \omega _{x} & 0 \end{array}\right), $$
(2)

where the indices ()s,g s have been omitted for readability. The accelerometer signal consists of a sensor acceleration component \(\user2{a}^{s}\) and a gravitational acceleration component \(\user2{g}^{s}: \user2{s}^{s} =\user2{a}^{s}-\user2{g}^{s}.\) The orientation \(\user2{R}^{g}_{s}\) is used to remove the gravitational acceleration \(\left(\user2{a}^{g} =\user2{R}^{g}_{s} \user2{s}^{s}+\user2{g}^{g}\right)\) [20], and the change of position in global coordinates is obtained by double integration of the sensor acceleration component \(\user2{a}^{g}.\) It should be noted that the orientation of the source with respect to the global \(\user2{R}^{g}_{c}\) is obtained in similar way using the acceleration \(\user2{s}^{c}\) measured by the accelerometer of the inertial sensor attached to the source, and its angular velocity \({\varvec{\omega}}^{c,g}_{c}\) measured by the gyroscope.

Fig. 1
figure 1

Overview of the measurement system used to estimate relative positions and orientations of an Inertial and Magnetic Measurement System (IMMS) with respect to the source. The source consists of three circular coils that are mounted orthogonally with respect to each other. An additional IMMS is mounted on the source to estimate its movement

Full size image
Fig. 2
figure 2

Relation between true sensor frame \(\Uppsi_{s},\) estimated sensor frame \(\hat{\Uppsi}_{s},\) magnetic source frame \(\Uppsi_{c},\) local coil frame \(\Uppsi_{cz}\) and global frame \(\Uppsi_{g}\) with a coil around the z axis. The relative position of the sensor with respect to the magnetic source is denoted by \(\user2{p}^{c}_{s},\) the relative position of the coil with respect to the magnetic source by \(\user2{p}^{c}_{cz},\) and the relative position of the sensor with respect to the coil by \(\user2{p}^{cz}_{s}\)

Full size image

The relative orientation of the sensor with respect to the magnetic source \(\user2{R}^{c}_{s}\) is obtained by

$$ \dot{\user2{R}}^{c}_{s} = {\user2{R}}^{c}_{s}\tilde{\varvec{\omega}}^{s,c}_{s}, $$
(3)

where the angular velocity is given by the difference between the angular velocities of the sensor and the source, both expressed in the same coordinate system

$$ {\varvec{\omega}}^{s,c}_{s} = {\varvec{\omega}}^{s,g}_{s} - \left({\user2{R}}^{g}_{s}\right)^{T} {\user2{R}}^{g}_{c} {\varvec{\omega}}^{c,g}_{c}. $$
(4)

Subsequently, the relative position \(\user2{p}^{c}_{s}\) is obtained by

$$ \begin{aligned} {\user2{a}}^{c}_{s} &= \left({\user2{R}}^{g}_{c}\right)^{T} \left({\user2{R}}^{g}_{s} {\user2{s}}^{s} - {\user2{R}}^{g}_{c} {\user2{s}}^{c} \right) = {\user2{R}}^{c}_{s} {\user2{s}}^{s} - {\user2{s}}^{c} \\ {\user2{v}}^{c}_{s} &= {\user2{v}}_{0} + \int\limits_{t_{0}}^{t}{\user2{a}}^{c}_{s}(\tau)\hbox{d}\tau \\ {\user2{p}}^{c}_{s} &= {\user2{p}}_{0} + \int\limits_{t_{0}}^{t}{\user2{v}}^{c}_{s}(\tau)\hbox{d}\tau. \end{aligned} $$
(5)

2.2 Filter structure

In this study, a complementary filter, commonly used for inertial navigation [2], is designed which operates on the errors in the state variables using a feedback structure (Fig. 3). The inertial measurements are not used purely as measurements, but rather as an input of the process model to form the reference trajectory against which the measurements of inclination and magnetic field are compared. The position and velocity change over time are extracted from measured inertial acceleration (Eq. 5) and the orientation is estimated from measured angular velocity (Eq. 3). It is the actual measurement minus the predicted measurement that forms the measurement error that is fed into the fusion filter. We chose to use an Extended Kalman Filter (EKF) that operates on the error states to fuse inertial with magnetic sensing. The error states for the filter are chosen as

$$ \delta {\user2{x}} =\left(\begin{array}{lllll} \delta {\user2{p}} & \delta {\user2{v}} & \delta {\varvec{\theta}} & \delta {\user2{a}} & \delta {\varvec{\omega}}\end{array}\right)^{T}, $$
(6)

where \({\delta}\user2{p}\) denotes the position error, \({\delta}\user2{v}\) the velocity error, \(\delta {\varvec{\theta}}\) the orientation error, \({\delta}\user2{a}\) the accelerometer bias, and \(\delta {\varvec{\omega}}\) the gyroscope bias. In general, the state space equations are given by

$$ \dot{\delta{\user2{x}}} = {\user2{f}}\left(\delta {\user2{x}}\right) + {\user2{w}}_{x} $$
(7)
$$ \delta {\user2{y}} = {\user2{h}}\left(\hat{\user2{x}}\right) - {\user2{y}} + {\user2{w}}_{y}, $$
(8)

where \(\user2{f}\left( {\delta}\user2{x} \right)\) denotes the function used to propagate the state vector \({\delta}\user2{x}\) in time, \(\user2{y}\) the measurement, and \(\user2{h}\left( \hat{\user2{x}} \right)\) the function which maps the estimated state vector \(\hat{\user2{x}}\) to the predicted measurement \(\hat{\user2{y}}.\) The additive process and measurement noises are denoted by \(\user2{w}_{x}\) and \(\user2{w}_{y}\) with covariances \(\user2{Q}_{x}\) and \(\user2{Q}_{y},\) respectively. The measurement model consists of a magnetic update \({\delta}\user2{y}_{m}\) and an inclination update \({\delta}\user2{y}_{i}.\) The inclination update is necessary, since an inaccurate estimation of inclination introduces integration drift caused by an incorrect removal of gravitational acceleration.

Fig. 3
figure 3

Structure of the Extended Kalman Filter (EKF). Measurements obtained from the Inertial and Magnetic Measurement System (IMMS) are used to estimate the position and orientation of the sensor with respect to the source \((\hat{\user2{p}}^{c}_{s}\) and \(\hat{\user2{R}}^{c}_{s})\) in an Inertial Navigation System (INS). A Kalman Filter (KF) uses an inclination update \({\delta}\user2{y}_{i}\) and a magnetic update \({\delta}\user2{y}_{m}\) to improve the estimation of \(\hat{\user2{p}}^{c}_{s}\) and \(\hat{\user2{R}}^{c}_{s}.\) Based on an estimation of the uncertainty of the position and orientation estimation \(\user2{P}^{-}_{k},\) the system decides if a magnetic update is needed. It also decides which coil needs to be actuated (coilnr) and which actuation current (I) should be used

Full size image

The recursive algorithm for a Kalman filter consists of a prediction step and an update step. Since the EKF resets the error states to zero immediately after a measurement, only the covariance is updated in the prediction step

$$ {\user2{P}}^{-}_{k} = {\user2{F}}_{k-1} {\user2{P}}_{k-1} {\user2{F}}_{k-1}^{T} + {\user2{Q}}_{x,k-1}, $$
(9)

where \(\user2{P}^{-}_{k}\) denotes the error covariance associated with the state \({\delta}\user2{x}\) at time instant k, and \(\user2{F}_{k}\) denotes the discrete process model that propagates the state in time. It should be noted that a minus superscript ()− denotes the a priori estimate. If a measurement becomes available, the state and covariance are updated resulting in the a posteriori estimate

$$ \begin{aligned} \delta \hat{{\user2{x}}}_{k} &= {\user2{K}}_{k} \delta {\user2{y}}_{k} \notag \\ {\user2{P}}_{k} &= \left({\user2{I}}-{\user2{K}}_{k}{\user2{H}}_{k}\right){\user2{P}}_ {k}^{-} \\ \delta {\user2{y}}_{k} &= \hat{{\user2{y}}}^{-}_{k} - {\user2{y}}_{k} \\ {\user2{K}}_{k} &= {\user2{P}}_{k}^{-} {\user2{H}}_{k}^{T} \left( {\user2{H}}_{k} {\user2{P}}_{k}^{-} {\user2{H}}_{k}^{T} + \user2{Q}_{y,k}\right)^{-1} \end{aligned} $$
(10)

with \(\user2{K}_{k}\) the Kalman gain, \(\user2{H}_{k}\) the linearized measurement model, and \(\user2{I}\) the identity matrix. Based on an estimation of the uncertainty associated with the position and orientation, represented on the diagonal of the covariance matrix \(\user2{P}_{k},\) the system decides when a magnetic update is necessary. In order to decide which coil needs to be actuated, the covariance update (Eq. 10) is calculated for each of the source coils that can be actuated. The system chooses to actuate only that coil with the highest contribution to the reduction of the uncertainty.

In Sect. 2.3, it is shown that the prediction of the state (Eq. 7) can be described by a linear equation using a single matrix \(\left(\user2{f}\left(\delta \user2{x}\right) =\user2{F}\delta \user2{x}\right).\) The nonlinear measurement models for the magnetic update \({\user2{h}}_{m}\left(\hat{\user2{x}}\right)\) and the inclination update \({\user2{h}}_{i}\left(\hat{\user2{x}}\right)\) as well as their linearized versions are described in Sect. 2.4.

2.3 Process model

The error equations that comprise the process model can be derived by writing the derivative of the state as a function of the state itself, according to Eq. 7. The error state \({\delta}\user2{x}\) denotes the difference between the true state \(\user2{x}\) and an estimation of the state \(\hat{\user2{x}}.\) For the position error, velocity error, accelerometer bias, and gyroscope bias, the relation is simply given by

$$ \hat{\user2{x}} = {\user2{x}} + \delta {\user2{x}}. $$
(11)

For the orientation error \(\delta {\varvec{\theta}},\) rotation matrices are used

$$ \hat{\user2{R}}^{c}_{s} = {\user2{R}}^{c}_{s} \delta {\user2{R}} \approx {\user2{R}}^{c}_{s} \left({\user2{I}} + \tilde{\delta {\varvec{\theta}}} \right) \Rightarrow \tilde{\delta {\varvec{\theta}}} \approx \left({\user2{R}}^{c}_{s}\right)^{T} \hat{\user2{R}}^{c}_{s} - {\user2{I}}. $$
(12)

The derivative of the orientation error is found by applying Eqs. 12, 3, and neglecting the product of errors

$$ \dot{\tilde{\delta {\varvec{\theta}}}} = \left(\dot{\user2{R}}^{c}_{s}\right)^{T} \hat{\user2{R}}^{c}_{s} + \left({\user2{R}}^{c}_{s}\right)^{T} \dot{\hat{\user2{R}}}^{c}_{s} \approx \tilde{\delta {\varvec{\omega}}} - \tilde{\varvec{\omega}}^{s,c}_{s} \tilde{\delta {\varvec{\theta}}} + \tilde{\delta {\varvec{\theta}}} \tilde{\varvec{\omega}}^{s,c}_{s}. $$
(13)

The vector representation of the orientation error \(\dot{\delta {\varvec{\theta}}}\) is found by inspecting the individual terms of the skew-symmetric matrix \(\dot{\tilde{\delta {\varvec{\theta}}}}\)

$$ \dot{\delta {\varvec{\theta}}} = \delta {\varvec{\omega}} - \tilde{\varvec{\omega}}^{s,c}_{s} \delta {\varvec{\theta}}. $$
(14)

The derivative of the velocity error is found by applying Eqs. 5, 11, 12 and neglecting the product of errors

$$ \dot{\delta{\user2{v}}} = \dot{\hat{\user2{v}}}^{c}_{s} - \dot{\user2{v}}^{c}_{s} = {\user2{R}}^{c}_{s} \delta {\user2{a}} - {\user2{R}}^{c}_{s} \tilde{{\user2{s}}^{s}} \delta {\varvec{\theta}}, $$
(15)

and the derivative of the position

$$ \dot{\delta {\user2{p}}}= \delta {\user2{v}}. $$
(16)

The accelerometer bias and gyroscope bias are modeled as first-order Markov processes

$$ \begin{aligned} \dot{\delta {\user2{a}}} &= - \beta_{a} \delta {\user2{a}} + \sqrt{2 \sigma^{2}_{a} \beta_{a} } {\user2{u}}\\ \dot{\delta {\varvec{\omega}}} &= - \beta_{\omega} \delta {\varvec{\omega}} + \sqrt{2 \sigma^{2}_{\omega} \beta_{\omega} }{\user2{u}}, \end{aligned} $$
(17)

where \(\user2{u}\) denotes unity white noise. Summarizing, the following state equations have been found

$$ \begin{aligned} \dot{\delta \user2{p}} &= \delta {\user2{v}}\\ \dot{\delta {\user2{v}}} &= {\user2{R}}^{c}_{s} \delta {\user2{a}} - {\user2{R}}^{c}_{s} \tilde{{\user2{s}}^{s}} \delta {\varvec{\theta}}\\ \dot{\delta {\varvec{\theta}}} &= \delta {\varvec{\omega}} - \tilde{{\varvec{\omega}}}^{s,c}_{s} \delta {\varvec{\theta}}\\ \dot{\delta {\user2{a}}} &= - \beta_{a} \delta {\user2{a}}\\ \dot{\delta {\varvec{\omega}}} & = - \beta_{\omega} \delta {\varvec{\omega}}. \end{aligned} $$
(18)

2.4 Measurement model

This section describes the measurement models that are used to update the estimation of the state. First, the model for the magnetic update \({\user2{h}}_{m}\left(\hat{\user2{x}}\right)\) and its linearized version \(\user2{H}_{m}\) are described, followed by the model for the inclination update \({\user2{h}}_{i}\left(\hat{\user2{x}}\right)\) and its linearized version \(\user2{H}_{i}.\)

2.4.1 Magnetic update

The magnetic measurement model is used to predict the field generated by the source coil at the location of the sensor based on an estimation of the position and orientation of the sensor. Figure 2 shows a configuration with a coil around the z axis, which means frames \(\Uppsi_{cz}\) and \(\Uppsi_{c}\) are aligned. The position of both frames with respect to each other is denoted by \(\user2{p}^{c}_{cz}.\) Applying the measurement model proposed in a previous study [18] to the configuration shown in Fig. 2 results in

$$ \hat{{\user2{y}}}_{m} = {\user2{h}}_{m}\left(\hat{\user2{x}}\right) = \hat{\user2{B}}^{s} = \left(\hat{{\user2{R}}}^{cz}_{s}\right)^{T} \hat{\user2{B}}^{cz}, $$
(19)

where \(\hat{\user2{B}}^{s}\) denotes the magnetic field measured by the sensor, and \(\hat{\user2{R}}^{cz}_{s} = \left({\user2{R}}^{c}_{cz}\right)^{T} \hat{\user2{R}}^{c}_{s}\) denotes the relative orientation of the sensor with respect to the coil. The magnetic field that is generated by the source coil at the location of the sensor expressed in coil coordinates is found by applying the Biot–Savart law to a circular wire loop (see Appendix), and expressing the results in cartesian coordinates

$$ \hat{\user2{B}}^{cz}(\hat{\user2{p}}^{cz}_{s}) = \frac{\mu_{0} N I} {2\pi\sqrt{\left(\sqrt{\hat{p}_{x}^2 + \hat{p}_{y}^2}+b\right)^{2} + \hat{p}_{z}^{2}}}\left( \begin{array}{l} \frac{\hat{p}_{x}\hat{p}_{z}} {\hat{p}_{x}^2 + \hat{p}_{y}^2}\left( -K(k) + \frac{b^2 + \hat{p}_{x}^2 + \hat{p}_{y}^2 + \hat{p}_{z}^2}{\left(\sqrt{\hat{p}_{x}^2 + \hat{p}_{y}^2}-b\right)^2 + \hat{p}_{z}^2}E(k)\right) \\ \frac{\hat{p}_{y}\hat{p}_{z}} {\hat{p}_{x}^2 + \hat{p}_{y}^2}\left( -K(k) + \frac{b^2 + \hat{p}_{x}^2 + \hat{p}_{y}^2 + \hat{p}_{z}^2}{\left(\sqrt{\hat{p}_{x}^2 + \hat{p}_{y}^2}-b\right)^2 + \hat{p}_{z}^2}E(k) \right) \\ K(k) + \frac{b^2 - \hat{p}_{x}^2 - \hat{p}_{y}^2 - \hat{p}_{z}^2} {\left(\sqrt{\hat{p}_{x}^2 + \hat{p}_{y}^2}-b\right)^2 + \hat{p}_{z}^2}E(k) \end{array}\right), $$
(20)

with μ0 the magnetic permeability of vacuum \((4\pi \cdot 10^{-7}\;\hbox{T}\;\hbox{m}^{2}/\hbox{A}),\) N the number of windings, I the current applied, and b the radius of the coil. It should be noted that the indices ()cz s have been omitted for readability. The relative position of the sensor with respect to the coil is given by \(\hat{\user2{p}}^{cz}_{s} = \left({\user2{R}}^{c}_{cz}\right)^{T} \left(\hat{\user2{p}}^{c}_{s} - {\user2{p}}^{c}_{cz}\right),\) and

$$ \begin{aligned} k &= \sqrt{\frac{4b\sqrt{\hat{p}_{x}^2 + \hat{p}_{y}^2}} {\left(\sqrt{\hat{p}_{x}^2 + \hat{p}_{y}^2}+b\right)^2 + \hat{p}_{z}^2}}\\ K(k) &= \int\limits^{\pi/2}_{0}\frac{1} {\sqrt{1-k^{2}\sin^{2}{\phi}}}\hbox{d}\phi \\ E(k) &= \int\limits^{\pi/2}_{0}{\sqrt{1-k^{2}\sin^{2}{\phi}}}\hbox{d}\phi. \end{aligned} $$
(21)

The linearized model is given by

$$ {\user2{H}}_{m} = \frac{\partial \hat{\user2{y}}_{m}}{\partial \delta {\user2{x}}} =\left( \begin{array}{llllll} \left(\hat{{\user2{R}}}^{cz}_{s}\right)^{T} & \frac{\partial \hat{\user2{B}}^{cz}}{\partial \delta {\user2{p}}} & {\user2{0}} & \tilde{\hat{\user2{B}}}^{s} & {\user2{0}} & {\user2{0}} \end{array}\right). $$
(22)

The partial derivatives of the estimated field with respect to the position error \(\left( \left(\hat{\user2{R}}^{cz}_{s}\right)^{T} \frac{\partial \hat{\user2{B}}^{cz}}{\partial \delta \user2{p}} \right)\) result in a rather lengthy and complicated expressions and are therefore not shown.

2.4.2 Inclination update

The inclination can be updated only during periods of low acceleration. Therefore, the measured acceleration is tested in advance for deviations from the gravitational acceleration. In case of significant deviations, the measurement noise is set to extremely high values such that the acceleration is not used to update the inclination during these intervals. The measurement equation for the inclination update is found by expressing the gravitational acceleration in sensor coordinates and adding the accelerometer bias

$$ \hat{\user2{y}}_{i}^{s} = {\user2{h}}_{i}\left(\hat{\user2{x}}\right) = \left(\hat{\user2{R}}^{g}_{s}\right)^{T} {\user2{g}}^{g} + \delta {\user2{a}} = \left(\hat{\user2{R}}^{c}_{s}\right)^{T} \left(\hat{\user2{R}}^{g}_{c}\right)^{T} {\user2{g}} ^{g} + \delta {\user2{a}}, $$
(23)

where \(\user2{g}^g\) denotes the gravitational acceleration in global coordinates.

The linearized model is given by

$$ {\user2{H}}_{i} = \frac{\partial \hat{\user2{y}}_{i}^{s}}{\partial \delta \user2{x}} =\left( \begin{array}{lllll} {\user2{0}} & {\user2{0}} & \tilde{\hat{\user2{g}}}^{s} & {\user2{I}} & {\user2{0}} \end{array}\right), $$
(24)

with \(\hat{\user2{g}}^{s}\) being the estimated inclination in sensor coordinates.

2.5 Experimental methods

The measurement system, as shown in Fig. 1, consisted of two IMMS modules (MTx with Xbus, Xsens Technologies B.V.), and three coils mounted orthogonal with respect to each other. Each coil was a circular coil with 50 windings and a radius of 0.055 m. The IMMS modules were sampled at 100 Hz. The maximum current that could be delivered by the driving electronics was 1.5 A. The driving electronics were controlled realtime by sending messages to the Xbus via a wireless Bluetooth connection. At certain time instants, the magnetic system generates magnetic pulses that were measured by the magnetometer of the IMMS. For each pulse, the mean and standard deviation were calculated which were both fed into the EKF. For validation, the relative position and orientation estimated by the fusion filter were compared to a reference optical position measurement system (Vicon, Oxford Metrics). Three markers with a diameter of 25 mm were attached to each IMMS. The reference data was obtained at a sample frequency of 100 Hz.

During the experiments, the system decided realtime when an actuation was required and which of the three coils should be actuated based on an estimation of the uncertainty associated with the position and orientation. The actuation current (1.3 A) and pulse duration (30 ms) were fixed. The complete experiment consisted of three parts. In the first experiment, the three source coils with IMMS were placed on a flat surface. An IMMS was moved around the source while varying the position and orientation. In the second experiment, the three source coils with IMMS were attached to the lower back of a subject, and the subject performed several movements with respect to the standing posture. With an IMMS attached to the back at the level of the first thoracic vertebra, the subject performed flexion/extension and rotation of the back. With an IMMS attached to the upper leg, the subject performed hip flexion/extension, and with an IMMS attached to the upper arm, the subject performed shoulder abduction/adduction. In the third experiment, the subject walked through the laboratory at a self-selected speed. The IMMS was attached to the upper leg and to the back at the level of the first thoracic vertebra.

3 Results

Figure 4 shows the relative position of the sensor with respect to the source for a representative trial of the first experiment, where the sensor was moved around the source coils which were not moving. The position error is defined as the difference between the estimated relative position by the ambulatory system and the reference system. Similarly, the orientation error is defined as the smallest angle about which the relative orientation of the sensor with respect to the source estimated by the ambulatory system has to be rotated to coincide with the relative orientation estimated by the reference system. The rms errors, averaged over 10 trials were calculated to be 0.028 ± 0.004 m (mean ± standard deviation) for the position and 3.1 ± 0.6 degrees for the orientation. A major part of the position error can be attributed to the peaks shown in the bottom figures of Fig. 4, which are caused by the experimental setup. During the peaks, a certain amount of data was queued in the input buffer causing the filter propagate slightly delayed data. If the filter then decides to perform a magnetic update, the input buffer must be emptied before the filter can process the response to the magnetic update causing the error to increase. The position error can be decreased by using the updated position after a magnetic actuation to reduce the drift between updates [20]. This resulted in a remaining rms position error of 0.022 ± 0.004 m averaged over 10 trials.

Fig. 4
figure 4

Relative position of the sensor with respect to the source for a representative trial of experiment 1. Top figures: x, y, and z coordinates estimated by the fusion filter (solid) and reference system (dashed). Bottom figures: position error

Full size image

As mentioned, the system decides realtime if it is required to actuate based on an estimate of the uncertainty associated with the relative position and orientation. Moreover, the system chooses which coil should be updated such that maximum information is obtained. An overview of the choices made by the filter for a trial, where the IMMS was moved around the source are shown in Fig. 5. Each dot represents a magnetic update and the coil that has been actuated is depicted by the value on the vertical axis. The right figure indicates that the choice of the coil to be actuated is indeed dependent on the location of the sensor.

Fig. 5
figure 5

The left figure shows the time instances and coil number of the coil that has been actuated. The right figure shows a top view of the relative position of the sensor with respect to the source. Each point indicates the position after an update with a source coil (x circle, y plus sign, and z cross). The origins of the source coils are indicated in the center of the figure (x circle, y plus sign, and z cross)

Full size image

During time intervals when the relative position and orientation of the sensor with respect to the source do not change, the time intervals between updates should increase. This is shown in Fig. 6, which shows the relative distance between source and sensor and the time instances of the magnetic updates. As can be seen, the time intervals between magnetic updates increase when the sensor is not moving, and decreases during movements. The bottom line indicates the uncertainty associated with the x position error \(\sigma_{\delta p_{x}},\) which remains bounded due to the magnetic updates.

Fig. 6
figure 6

Decreasing update rate when sensor is not moving. The relative distance between sensor and source is indicated by the solid line for the fusion filter and the dashed line for the reference system. The time instances of the magnetic updates are indicated by the dots. The bottom line indicates the uncertainty associated with the x position error \(\sigma_{\delta p_{x}}\)

Full size image

An example of the relative position of the sensor with respect to the source for a shoulder abduction/adduction trial with an IMMS attached to the upper arm is shown in Fig. 7. The solid line indicates the estimation by the fusion filter, the dashed line indicates the estimation by the reference system. The rms position error was 0.025 ± 0.012 m, the rms orientation error was 2.8 ± 0.7 degrees averaged over five trials.

Fig. 7
figure 7

Relative position of the sensor with respect to the source for shoulder abduction/adduction with an IMMS attached to the upper arm, estimated by the fusion filter (solid) and reference system (dashed)

Full size image

The bottom figure of Fig. 8 shows the relative position of the sensor with respect to the source for a representative walking trial with the IMMS attached to the back at the level of the first thoracic vertebra. Despite the movement of source and sensor as indicated by the top figure of Fig. 8, the filter is able to provide a stable estimate of the relative position and orientation. The rms errors, averaged over six trials, were calculated to be 0.026 ± 0.004 m for the position and 3.6 ± 0.4 degrees for the orientation. A complete overview of the position and orientation errors for all experiments that were performed is shown in Table 1.

Fig. 8
figure 8

Representative trial of a subject during walking with the source attached to the lower back and an IMMS at the level of the first thoracic vertebra. Top figure: distance traveled by the subject estimated by the reference system. Bottom figure: relative position of the sensor with respect to the source estimated by the fusion filter (solid) and reference system (dashed)

Full size image
Table 1 Mean and standard deviation (SD) of rms position and orientation error for all movements that were performed during the experiments
Full size table

4 Discussion

The present study leads to the conclusion that the proposed adaptive filter allows accurate ambulatory tracking of relative positions and orientations on the human body. The change of position and orientation is estimated using inertial sensing. A fusion filter (EKF) estimates the uncertainty associated with the position and orientation and a magnetic coil is actuated only if the uncertainty associated with the position or orientation exceeds a predefined threshold. Moreover, the coil delivering most information is actuated only. The actuation current and pulse duration were fixed during the experiments. In future experiments, these parameters can also be tuned adaptively such that maximal accuracy is achieved at minimal energy consumption.

The results indicate that the filter is able to provide a stable estimate of relative position and orientation for several types of movements. Although a previous study [14] showed higher accuracy for the relative position (position error approximately five times smaller), and comparable accuracy for the relative orientation, significant improvements were achieved for the update rate (average time between magnetic actuation of 547 ms using single coil actuation instead of 600 ms using three coils actuation), and pulse width (30 ms instead of 60 ms). Moreover, the orientation of source and sensor reported in [14] were obtained from a separate Kalman Filter, whereas this study estimates the orientation of source and sensor in a single filter running realtime. In order to achieve higher accuracy, it is suggested to improve the experimental setup such that the cause of the error peaks reported in Sect. 3 is removed. An improved setup should also allow multiple coils to be actuated simultaneously during rapid movement periods, which is not possible using the current setup. Another suggestion to improve the accuracy is to apply a Kalman smoothing algorithm in an off-line analysis that propagates the filter backward in time [9].

The presence of nearby ferromagnetic materials can influence the accuracy of the filter negatively. In our previous study [18], it was shown that a ferromagnetic material influences the measurements only if it was located near the sensor or between source and sensor. It is therefore important to remove any ferromagnetic objects from the body during the measurements. During periods of no actuation, the magnetometers can be used to detect ferromagnetic disturbances. In an undisturbed environment, the magnetometer will measure a homogeneous earth magnetic field. In case of a nearby ferromagnetic object, the magnetometer can be used to detect the inhomogeneity of the magnetic field caused by the ferromagnetic object. By increasing the measurement noise of the magnetic update during these periods, the filter relies less on the disturbed measurement.

Despite the promising results, several aspects of the proposed system can be improved. An improvement could be to reuse the energy of the magnetic field for future actuation. The hardware used for the experiments is designed such that it dissipates all energy needed to build up the magnetic field. During the pull down phase of the current pulse, the change of magnetic field can be used to induce a current in another source coil. The electrical energy can be stored and used to build up the field for the next actuation.

Another improvement would be to use an adaptive instead of a fixed threshold to decide if the magnetic system should actuate, by taking sensor location into account. If the sensor is located near the source, a relatively small amount of energy is required for an update. If the sensor is located far from the source, much energy is required for an update while the reduction of the uncertainty will be small since the signal to noise ratio will be low. Since the magnetic field decreases with the third order of distance and the magnetic field is proportional to the current applied, a magnetic field with equal magnitude at an increased distance yields an electric current increase of the third order. This means that by adapting the uncertainty threshold and actuation parameters based on the relative distance between source and sensor such that the system increases the update rate for sensor locations near the source, the overall energy consumption can be decreased.

It should be noted that other effective solutions exist to obtain an estimate of relative orientation of sensors with respect to each other. Accelerometers provide an estimate of inclination during periods of low acceleration (Sect. 2.4.2). Similarly, magnetometers provide an estimate of heading using the earth magnetic field [13, 16, 23]. The heading accuracy depends on the amount of earth magnetic field disturbance [4]. Moreover, the earth magnetic field does hardly provide heading information if the direction of the magnetic field is close to the vertical, which may occur in a disturbed environment. A major advantage of the system proposed is that it does not need the earth magnetic field for an accurate estimation of relative position and orientation.

Recently, we proposed an instrumented shoe that can be used to assess ankle and foot dynamics [20], and center of mass movement during walking [19] in an ambulatory environment. An important aspect missing in the instrumented shoe principle is the estimation of relative positions, for example, the distance between the feet or the distance between the center of mass and the feet. Although these distances cannot be covered using the setup used in the present study, the coil dimensions can be optimized to be suitable. Concluding, the ambulatory tracking system is expected to provide a valuable contribution to human movement analysis, as it allows relative positions and orientations to be estimated using a wearable system.