1. Introduction
In the last decade advances in microelectromechanical sensors (MEMS) have propelled biomechanics applications based on inertial measurement units (IMUs) forward [
1,
2,
3]. Many of these applications are based on estimated IMU orientation. Orientation estimates are usually obtained through fusion of the triaxial signals of accelerometers and gyroscopes [
4].
Despite the many advances, low cost MEMS used in biomedical applications still experience bias, scale factors and other random errors [
4,
5]. These errors are due to misalignment of the axes in manufacturing of the sensor or during assembly of the IMU, and due to self-heating effects [
5,
6]. Adequate sensor calibration and fusion methods can reduce or compensate the effect of these errors. In calibration a relation is established between the digital sensor measurements and the physical quantity being measured. Calibration methods differ mainly in the instrumentation used and the complexity of the model used to represent the sensor (see [
5] for an overview).
The most common fusion method for IMU orientation is Kalman Filtering (KF) [
4]. In KF, the primary information is obtained through integration of the gyroscope signals. Errors however accumulate quickly when integrated, resulting in a drift that exponentially deteriorates the orientation estimate. In many KFs accelerometers therefore fulfill the role of inclinometers, providing the initial orientation and acting as secondary sensors to correct the gyroscope based estimate [
4]. Accelerometer measurements consist of accelerations due to both gravity and body motion. Gravitational acceleration must dominate the accelerometer measurement in order for it to be used as an inclinometer. However, accelerations due to human motion are time variable, segment dependent, and can be high under dynamic conditions. It is therefore important to either seperate both sources of acceleration [
7,
8], or accurately establish a criterion to determine the absence of acceleration due to human motion [
7,
9]. For clarity we have chosen to only compare KFs using a criterion approach to determine accelerometer reliability.
KFs can also be grouped based on how they model the state-space. Local KFs only incorporate information from the individual IMU in their model, and can thus only use information from that IMU itself to improve its orientation estimate. In combination with criterion based accelerometer screening, this can lead to an absence of updates under dynamic conditions. This is especially true in segments with little static periods, such as the trunk. The authors have therefore previously presented global or cooperative KF models [
10,
11]. The cooperative model assumes a relation between multiple sensors and uses the information present in the collective to update the individual sensors. In [
10] information was shared between the most reliable IMU, in that time frame, and the potentiometers embedded in the exoskeleton, using a Markovian Jump Linear System.
Inertial sensors attached to the trunk are commonly used to obtain orientation of the exoskeleton. In [
11], information was shared between all the IMUs located on the lower limb. The latter approach is independent of the exoskeleton data and made use of all IMU data instead of only the most reliable IMU. It is unclear which of these cooperative methods is more susceptible to noise in the IMU measurements.
Most biomechanics applications and studies are performed with commercial off-the-shelf IMUs (e.g., Xsens, Technaid). These high-end IMUs are calibrated by the manufacturer and provide orientation based on an embedded fusion algorithm. High-end IMUs can cost over 1000$ per unit and rarely provide real-time access (e.g., through CAN-bus) to the raw sensor data. The high price and closed sensor fusion and communication algorithms make these high-end IMUs appear less attractive to be used embedded in a device or platform, e.g., an exoskeleton. Recently, many sensors and modules have become available that can present an alternative to high-end IMUs for certain applications [
12,
13]. These sensors however, often are not calibrated when mounted or were subjected to a minimal or limited calibration protocol. Many applications and studies only take into account the physical data, assuming the transformation from digital to physical format to be correct [
5]. This calibration step is not trivial, in particular for gyroscopes [
5]. Despite a significant body of literature dedicated to both sensor calibration and sensor fusion, currently no studies have investigated their interaction. It is unknown whether limitations from simplified calibrations protocols can be aleviated by cooperative filtering, and whether the gains of more complex calibrations affect cooperative and local KFs equally.
In this paper we compare the influence of various calibration methods on different nominal Kalman Filtering (KF) algorithms. The objective is to demonstrate the gains that can be achieved by using different calibration and/or filtering methods, and to provide recommendations on optimal combinations for different situations. In order to do so, we validate both methods in the calibration with an exoskeleton. This allows us to use the sensors embedded in the exoskeleton for the Markovian Jump Linear System Kalman Filter. To facilitate the comparison between the KFs used in this study, all KFs used apply the same criterion approach to determine the reliability of the accelerometer data.
2. Experimental Section
Reliability of absolute angles estimation is important in both gait analysis and exoskeleton control. In these applications, absolute angles represent the orientation of a segment with respect to the global reference, in inertial sensors the earth frame. In
Figure 1 the absolute angles of the trunk (
), thigh (
), shank (
), and foot (
) are shown. Absolute segment angles are estimated using:
with
being the estimated absolute angle based on the integrated gyroscope data, and
, for
, the correction estimated by the KF based on the accelerometer data.
In this study, two calibration methods are applied in combination with three nominal KFs. One discrete local KF and two cooperative KFs are compared. All KFs are based on the same threshold criterion Ψ to determine realiability of the accelerometer data [
14,
15]. When the criterion is not met, only the prediction steps of the KF are performed. Updates are thus only performed when the accelerometer data is deemed reliable according to:
with
being the norm of the measured acceleration,
the magnitude of the gravitational acceleration, and
ζ being an arbtirarily small value between 0 and 1, depending on the IMU used. The Local KF’s model is only based on the individual IMU, whereas the cooperative KFs take advantage of the information available in the collective of sensors. The cooperative KFs used in this paper are a Markovian KF [
10] and a Matricial KF [
11]. The Markovian KF combines inertial sensor data with data from the potentiometers of the exoskeleton. In the Matricial KF only data from the IMUs is combined. We will first describe the protocol and instrumentation used. Subsequently, we describe the two calibration protocols,
simplified and
complex calibration, followed by a brief explanation of the KFs.
2.1. Protocol and Instrumentation
Data from a single healthy subject (33 years old male, 1.75 m, 82 kg) are presented. The subject had ample experience in walking with the exoskeleton and walked for 60 s on a treadmill at 1.5 km/h. The exoskeleton assisted his gait following a prescribed trajectory using a position based control algorithm. All systems used were synchronized using an electronic trigger. The first and last 10 s of the trial were discarded. The remaining 40 s were seperated into 20 s for optimisation of the filters, and 20 s for validation. Filter parameters were optimised using a genetic algorithm presented in [
16] (see online supplementary material for the parameter values). The point at which a KF is initialized is known to influence the convergence and thus performance of the KF. Therefore, six different starting points were chosen in both optimisation and validation phases (
Figure 4). This resulted in 36 combinations for each filter.
Kinematic data was acquired from three independent systems: an optoelectronic system (BTS Bioengineering, Italy), potentiometers embedded in the H2-exoskeleton (Technaid S.L., Madrid, Spain), and inertial sensors (Technaid S.L., Spain). Both the passive markers of the optoelectronic system and the inertial sensors were mounted on rigid boards that in turn were attached to the exoskeleton (
Figure 1). Care was taken to align all systems in the sagittal plane prior to data-collection. An additional alignment in the sagittal plane was performed post-data collection through a custom Matlab script. The kinematic data from the optical system was fitted with a spline. The orientation of the segment was subsequently obtained using the Gram-Schmidt method. The optical system global vertical axis is assumed to be aligned with gravity. A residual matrix was calculated to align the optical data to the IMU data [
17].
The reference measurement was taken from the optoelectronic system. Clusters of three markers were used to compute the segment orientation using the Gram-Schmidt method of consecutive cross-products. The clusters used in the data-analysis are highlighted in red in
Figure 1. Clusters were selected based on minimal marker occlusion during the trial. Optoelectronic data was acquired at 100 Hz and downsampled to 50 Hz using cubic spline interpolation.
Four IMUs were used, attached to the right leg and pelvis. The foot sensor was placed on top of the foot, aligned with both the sagittal plane of the leg and the longitudinal axis of the foot. The IMUs placed on the shank, thigh and pelvis were placed directly on the exoskeleton in the sagittal plane (
Figure 1). To ensure alignment, all sensors were rotated to be perfectly aligned in the sagittal plane via a static calibration. IMU data was acquired in raw digital format at 50 Hz. This data was then processed offline using custom Matlab scripts. First, the data was transformed from digital to physical format based on above mentioned calibration files. This transformation was repeated once for each of the four different calibration files. The physical data was subsequently used in the different KFs. The parameters of each of the filters were optimized based on the first 20 s of data using a genetic algorithm as presented in [
10,
16].
The H2-exoskeleton is a bilateral exoskeleton with six actuated degrees of freedom in the sagittal plane [
18]. The H2 is designed for gait rehabilitation of adults with a stature between 1.50 m and 1.90 m, and a bodyweight below 100 kg. Users of the H2 must have trunk stability while maintaining standing balance. The maximum walking speed the H2 can attain is 2km/h; the trial was performed at 1.5 km/h. The embedded potentiometers have a reported linearity of 0.25% [
19] and measure the relative angle between the segments of the exoskeleton. These relative angles correspond to the sagittal plane angles of the hip, knee, and ankle. The H2-potentiometer data was acquired at 1 kHz. This data was then downsampled to 50 Hz using cubic spline interpolation, and processed in Matlab.
2.2. Data Analysis
Performance of the KF in combination with the various calibration routines was assessed based on mean error (ME) and root mean square error (RMSE) metrics:
The statistical analysis aimed at assessing the effect of the calibration and filters on the various segments. A three-way mixed ANOVA was applied, with calibration and filter as within-subjects factors. The depedent measurement being the RMSE. Tukey-LSD post-hoc testing was perfomed where needed, with Bonferroni correction for multiple comparisons. The statistical analysis was performed using IBM SPSS Statistics software package (IBM SPSS Statistics 23, SPSS IBM, New York, NY, USA).
2.3. Calibration
The TechMCS IMUs (Technaid SL, Madrid, Spain) used in this study contain three types of sensors: accelerometers, gyroscopes, and magnetometers. Two calibration routines were applied to the digital data acquired from these sensors. Each calibration routine was repeated two times. All models applied in calibration are linear sensor models. The instrumentation and protocol used differs between each of the calibrations. The KFs applied in this study only rely on the data from the accelerometers and gyroscopes. Since the calibration of the former is relatively straightforward, the diference between both methods is predominantly related to the gyroscope calibration.
The setup of the first calibration method, the
simplified calibration, is minimal. It should be possible to replicate this simplified calibration in most laboratories. This calibration could thus easily be applied to IMUs that have been mounted or purchased without prior calibration. The accelerometers were calibrated placing the sensor housing level for each axis, so that it aligns with gravity (for both positive and negative directions: X, -X, Y, -Y, Z, -Z). The gyroscopes were calibrated using 90 degree angle executions. We used a custom mechanical platform to generate the 90 degree angles for each axis and direction (
Figure 2). By spinning the wheel, consecutive 90 degree angle executions were performed by the IMU. The gyroscope data was corrected by least-squares fitting of the integrated angular velocity signal to the 90 degree reference. Temperature effects are taken into account by calibration in a heat chamber, where gradual temperature changes between 10 and 50 degrees Celsius were generated. A linear temperature model was applied, as suggested in [
6].
In the
complex calibration method, the custom manual platform used for gyroscope calibration in the simplified calibrations is replaced by a high precision turntable (AC11205, Acutronic, Hirzel, Switzerland) with a rate stability of 0.0001% over 360 degrees (
Figure 3). This setup is closer to what might be used in high-end calibrations. In previous work we demonstrated that there is a velocity dependent component to the error in gyroscope signals [
20]. In the complex calibration, each IMU was therefore rotated about each of its axis at various angular velocities. The velocities applied were positive and negative 100, 200, 300, 400, 500 degrees/second. The measured angular velocity was fitted to the reference angular velocity provided by the turntable using a least-squares approach.
2.4. Filters
A. Local KF
The discrete local KF is the most common KF implementation to estimate absolute orientation. The model used in this KF is considered local since it is only based on the individual IMU, no information from other sensors is used. The state space model and output equation of the local KF are:
with
and
being, the orientation estimation and offset errors of the gyroscope;
and
the white noise Gaussian of the gyroscope and the gyroscope bias, and
the Markov process correlation time [
10,
16].
Based on the model and output equation, the KF can estimate the state, represented by
. In Algorithm 1 the discrete local KF is presented with following matrices [
21,
22]:
,
,
.
With ζ being an arbtirarily small value between 0 and 1.
In cooperative KFs, Equations (5) and (6) are written in matrix form:
The state is then given by:
where
;
, are the errors between the absolute angles (
) and the angle estimates calculated by the gyroscopes (
). Information from the collective of sensors is thus used to improve the absolute angle estimate of each segment.
Considering the gyroscope and gyroscope bias models (see Equation (5)), matrix
is represented by:
The time-invariant weighting matrix
Q is given by:
The difference between both types of models is most evident in the
Update equations. Consider a simplified scenario with a body consisting of two segments A and B. In the local scenario the correction terms would be:
As can be seen in Equation (
10), the correction term
is thus only based on data from the sensor attached to segment A. Whereas in the cooperative scenario the correction terms
and
are given by a weighting of both segments:
B. Cooperative Matricial KF
The Matricial KF fuses the information from all IMUs [
11]. To ensure accelerometer reliability, whilst maximizing the number of sensors used in the update, a minimum of two IMUs have to provide reliable accelerometer data before an update is performed. This is in agreement with the findings from [
11]. Using four IMUs the matrix
is then defined as:
which produces an output vector
The time-invariant matrix
R is defined as
C. Cooperative Markovian KF
The Markovian KF only uses the information from the most reliable IMU, and combines this with the information from the potentiometers from the H2-exoskeleton [
10].
The most reliable IMU is the IMU that passes the criterion displaying the lowest acceleration. The criterion for reliable accelerometer measurement and the Markovian state are defined as:
where
describe the lower limb or exoskeleton segments,
is the current Markovian state and
is an index that describes the reliability of the accelerometer used at the Markovian state,
. Changing from using one IMU to another IMU only depends on the current timeframe and is considered a jump. This filter is therefore also referred to as a Markov Jump Linear System.
3. Results and Discussion
Table 1 reports the RMSE obtained between each combination of filter and calibration applied in this study. A consistent trend to a decrease in RMSE can be seen between the simplified and more complex calibrations. The simplified calibrations furthermore appear to have a higher variability. Overall the thigh and shank absolute angles were quite robust to changes in calibration quality and show good accuracy. The trunk segment was the most critical, with RMSE as high as 4 degrees. The Markovian KF appeared to perform the best. The bottom row shows the percentage of total trial duration the accelerometers are used by the respective filters. The mean percentage reported for the Local KF, is the mean across all four segments. This mean however hides a large variability. For example, under calibration 1 the shank is updated 32% of the trial duration but the trunk only 7% of the trial duration. A more detailed interpretation of the differences and trends observed in
Table 1 is provided through the statistical analysis of the data.
In
Figure 4 the spread of the updates across the trial is indicated by coloured circles. Each coloured circle represents a KF update, with the colour indicating the applied calibration file. The pink, green, purple and yellow dots respectively stand for the data obtained applying calibration 1, 2, 3, and 4. Calibrations 1 and 2 refer to the simplified calibration, calibration 3 and 4 refer to the complex calibration. The light grey lines represent the absolute angle obtained from the accelerometers, whereas the dark grey line represents the absolute angle estimate based on the gyroscope data. An integration drift can be observed in the gyroscope based estimate. Both accelerometer and gyroscope based estimates form the input data to the KF algorithms. The different filter initiation points, and their respective intervals, are marked with horizontal arrows in
Figure 4A. The panels represent from top to bottom, the local KF, the Markovian KF, and the Matricial KF. All data shown in
Figure 4 relates to the trunk. In
Figure 4A, it is clear that, even at a relatively slow walking pace of 1.5 km/h, applying the Local KF results in few updates for the trunk. The Matricial KF (
Figure 4C) has a higher number of updates, but they appear to be concentrated in the stance phase. Few updates occur during swing. The Markovian KF (
Figure 4B) has a lower total number of updates than the Matricial KF (see also
Table 1), but they are spread more homogenously across the trial.
Figure 5 shows the performance of the three filters for each of the calibrations, for the trunk segment. The reference absolute trunk angle is represented by the discontinous black line. Absolute angle based only on accelerometer data (light grey), or gyroscope data (dark grey) is also shown to demonstrate the necesity and efficacy of sensor fusion.The panels again represent the Local KF (
Figure 5A), the Markovian KF (
Figure 5B), and the Matricial KF (
Figure 5C). The pink, green, purple and yellow lines respectively stand for the data obtained applying calibration 1, 2, 3, and 4. Calibrations 1 and 2 are repetitions of the simplified calibration routine. Calibrations 3 and 4 are repetitions of the more complex calibration scheme. The coloured lines represent the average data of the 6 validation datasets for that segment. The thicker the coloured line, the higher the variance present in the estimate. The Local KF appears to perform the poorest. Large errors with respect to the reference can be observed in the righthand side of
Figure 5A. The performance of the Markovian KF appears to contain more variance, while the divergence from the reference is more present in the pink and yellow line of the Matricial KF.
A three-way mixed ANOVA was performed to assess the effect of the calibrations on the filters and their interactions. No significant differences between calibrations 1 and 2, and calibrations 3 and 4 were observed. The analysis was therefore performed pooling the simplified calibrations (1 and 2) and the complex calibrations (3 and 4). The main effects of calibration, filter, and segment were all significant at the 5% level. A post-hoc testing revealed that the simplified calibration performed significantly poorer than the more complex calibration, irrespective of the filter applied or the segment under investigation. There was a main effect of segment (F(3, 105) = 151.7 ,
p = 0.000), as shown in
Figure 6. This effect tells us that if we ignore the type of filter and the calibration that were applied, the RMSE of the absolute angle estimates of the segments were significantly different. This was expected, since the trunk is a segment with fewer quasi-static moments in which the accelerometer data can be used to perform the update. The foot segment on the other hand is characterised by higher peaks related to initial and final contact. Post-hoc analysis indicated that the trunk performed significantly worse compared to the other segments. The foot was also significantly different from the shank and thigh. No significant difference was found between the shank and the thigh.
Figure 7 reports the main effect of filter. The Mauchly’s test of sphericity was significant (
p = 0.000,
ϵ = 0.793), the Greenhouse-Geisser correction was therefore applied to the data (F(1.58, 55.3) = 4.37,
p = 0.021). Post-hoc testing revealed a significant difference only between the Markovian KF and the Matricial KF. This effect tells us that if we ignore the segment and the calibration, the Markovian filter performed significantly better then the Matricial KF.
The interaction between the filter and the calibration is reported in
Figure 8. This interacion is significant (F(2,70) = 4.40,
p = 0.00). This effect tells us that the influence of applying different filters was different in the simplified calibration compared to the complex calibration. For each filter there is a significant improvement when applying the complex calibration compared to the simplified calibration. The Markovian KF outperforms the other filters when the simplified calibration is applied. This difference is however not significant. When the complex calibration is applied drastic improvements are observed, in particular for the Matricial KF. The Matricial KF performs significantly better than the Local KF. The significant main effect of filter between the Markovian and Matricial KF (see
Figure 7) appears to be due to the poor performance of the Matricial filter under the simplified calibration.
In
Figure 9 the interaction between the calibration and filter is shown for the trunk segment. Isolating the trunk, all filters improve significantly when changing from a simplified calibration to a complex calibration. The Matricial KF in particular improves from being the worst to the best performing filter. It appears that the error present in the individual sensors when using the simplified calibration protocol are augmented by this filter. When the data is of higher quality, the Matricial KF is able to benefit more than the Markovian KF. No significant difference between the Matricial and Markovian KF is observed when using the complex calibration.
4. Conclusions
This paper presented a comparison between local and cooperative Kalman Filters to estimate the absolute segment angle. The filter performance was assessed under two sensor calibration conditions. A simplified calibration, that can be replicated in most labs; and a complex calibration routine, with a setup similar to that used by high-end IMU manufacturers. Walking data from a healthy subject wearing an exoskeleton was collected. The subject walked for one minute at constant pace of 1.5 km/h on a treadmill. The hypothesis was that cooperative KFs would perform better than a local KF, in particular as the quality of the calibration improved. Significant effects of filter and calibration, as well as a significant interaction effect between filter and calibration were observed. Since there was also a significant effect for segment, the analysis of the paper describes both the general effect and the effect on the trunk. The trunk was identified from literature as being the most sensitive to changing the applied KF. This was confirmed by the results in this paper. The Markovian KF perfomed better than the Local KF, irrespective of the IMU calibration. If data from external sensors, such as potentiometers, are available incorporating them into the filter can result in an improved performance, especially if the calibration is less ideal or unknown. Our results further suggest that IMU calibration significantly affects performance, irrespective of the applied filter. The results indicate that where possible, IMUs with a high-end calibration should be used.
The Matricial KF was the most sensitive to IMU calibration, displaying the biggest performance gain when changing from the simplified to the complex calibrations. When using sensors that have been subjected to a high quality calibration, incorporating the potentiometer data thus does not result in a significant gain compared to using a cooperative filter, such as the Matricial KF. In future work we will add in-field calibrations to the comparison, as well as other less repetitive, tasks such as reaching and grasping and object manipulation. We will furthermore also evaluate cooperative KF with a more explicit physical model declaring the relation between the IMUs.