Identification of Noise Covariance Matrices to Improve Orientation Estimation by Kalman Filter

Magneto-inertial measurement units (MIMUs) are a promising way to perform human motion analysis outside the laboratory. To do so, in the literature, orientation provided by an MIMU is used to deduce body segment orientation. This is generally achieved by means of a Kalman filter that fuses acceleration, angular velocity, and magnetic field measures. A critical point when implementing a Kalman filter is the initialization of the covariance matrices that characterize mismodelling and input error from noisy sensors. The present study proposes a methodology to identify the initial values of these covariance matrices that optimize orientation estimation in the context of human motion analysis. The approach used was to apply motion to the sensor manually, and to compare the orientation obtained via the Kalman filter to a measurement from an optoelectronic system acting as a reference. Testing different sets of values for each parameter of the covariance matrices, and comparing each MIMU measurement with the reference measurement, enabled identification of the most effective values. Moreover, with these optimized initial covariance matrices, the orientation estimation was greatly improved. The method, as presented here, provides a unique solution to the problem of identifying the optimal covariance matrices values for Kalman filtering. However, the methodology should be improved in order to reduce the duration of the whole process.


Introduction
In the fields of medicine and biomechanics, human motion analysis is often used to extract quantitative and objective data to characterize a movement or a function of the neuro-musculo-skeletal system. During human motion analysis, the aim is frequently to define a body segment rotation with respect to its proximal segment [1]. To perform this analysis, optoelectronic systems that measure the 3D coordinates of reflective markers remain the "gold standard" [1,2]. However, such systems suffer from four main drawbacks: (a) even if such systems could be movable, because cameras have to be handled with caution, most systems are fixed inside the laboratory walls; (b) measurement volume is restricted to the camera's field of vision; (c) sessions with multiple subjects require numerous cameras to counteract occlusion by the marker; (d) calibration of the whole system can be time-consuming and must be redone if a camera has been moved.
In this context, magneto-inertial measurement units (MIMUs) are a promising, emerging way to perform human motion analysis outside the laboratory, because they are wearable, light, and wireless [3,4]. These MIMUs typically include a combination of 3D accelerometers, gyroscopes, and magnetometers. However, MIMUs do not measure orientation directly. Rather, they measure its linear acceleration, angular velocity, and the magnetic field surrounding the sensor. Angular velocity can be numerically integrated to obtain rotation, but MIMUs designed for human motion analysis are based on low-cost microelectromechanical technology that also suffers from substantial noise and bias over the measurement [5]. This noise and bias particularly affects the angular velocity measure, [6] such that the estimation of the orientation by numerical integration leads to a well-known error: orientation drift [7]. To improve the orientation calculation, the sensor's data are then combined [8], the linear acceleration and the magnetic field being used to correct the initial estimate obtained by numerical integration of the angular velocity [7]. Indeed, in the absence of movement, the linear acceleration measured by the sensor should be gravity, which is known to be vertical. This information can thus be exploited to define the inclination of the sensor. The measured magnetic field can also be used to express the orientation of the sensor relative to the Earth's magnetic field. However, in practice such ideal conditions are never obtained. Since the MIMUs are used for the purpose of characterizing human movements, accelerometers measure external acceleration induced by body segment movements in addition to the acceleration induced by gravity. In addition, magnetometers are sensitive to the ferromagnetic components located near the MIMU, which is a recurrent issue.
The classic approach to processing all of these noisy and perturbed measurements to estimate MIMU orientation is to fuse them inside a Kalman filter. In fact, Kalman filtering is largely adopted to perform data fusion from noisy measurements [8,9]. Complementary filters have been more recently proposed as an alternative to Kalman filtering [10,11]. Some studies have compared the performance of such filters with that of Kalman filters during human movement [12,13]. Since the performances were relatively similar, further studies must be performed to assess the advantages and disadvantages of these two types of filtering, by measuring various human movements in terms of duration and movement acceleration. The principal concept of Kalman filtering is to combine different sensor measurements with a dynamic model that translates the time-behavior of the state to be estimated. Despite its extensive use in MIMU processing, there is no consensus as to the structure of the Kalman filter and its settings [8]. In particular, a critical point when implementing a Kalman filter is the definition initialization of the covariance matrices that characterize mismodelling and input error from noisy sensors [14]. These covariance matrices are crucial since their poor definition can lead to inaccurate estimation as well as divergence and stability problems [15]. When applied to the orientation estimation by MIMU data, covariance matrices represent measurement noise from gyroscopes, accelerometers, and magnetometers. However, they also represent the external perturbations of acceleration and the magnetic field since, as previously mentioned, both the external acceleration experienced by the sensor, and magnetic disturbances, tend to cause errors in the estimated orientation. The critical issues of the initialization of covariance matrices are little documented [16,17] and, surprisingly, are never considered in human motion captured by MIMUs.
The aim of the present study is to propose a methodology which will identify the values of these covariance matrices that optimize orientation estimation in the context of human motion analysis.

Materials and Methods
To identify the optimal values of the covariance matrices, we proposed to explore the space of sets of values for each parameter of the covariance matrices, and to compare Kalman filter results with a reference measurement. For these purposes, we compared the orientation of a rigid object measured simultaneously by an optoelectronic system and by an MIMU. The optoelectronic system was considered in this case to be the most suitable system, due to its faculty to measure a wide range of movements in terms of acceleration and angular velocity level.

Experimental Procedure
The optoelectronic measurement was provided by a system of 20 cameras (Vicon Motion Systems Ltd., Oxford, UK) with a sampling frequency of 250 Hz. The object used was a set square, seen in  Figure 1, the usual function of which is to calibrate the optoelectronic system. This object was particularly suitable for this study because it is pre-equipped with five reflective markers that can form a coordinate system based on the geometry of the set square. Additionally, the handle makes it easy to apply hand movements, thus limiting marker occultation.
The MIMU (Opal sensors, APDM, Portland, OR, USA) was aligned with the set square by means of a custom adapter. Sensor data (triaxial accelerometers, gyroscopes, and magnetometers) were collected at 128 Hz.
Time synchronization between both systems was ensured by different procedures. During the experiment, an external trigger executed by the launch of the MIMU system measure triggered the optoelectronic system measure. However, after resampling the measurement from the optoelectronic system by cubic interpolation at 128 Hz to be consistent with the inertial measurement, post-processing showed a time-shift between the data. Therefore, another procedure, depicted in Annex A, was applied to time-synchronize the data. seen in Figure 1, the usual function of which is to calibrate the optoelectronic system. This object was particularly suitable for this study because it is pre-equipped with five reflective markers that can form a coordinate system based on the geometry of the set square. Additionally, the handle makes it easy to apply hand movements, thus limiting marker occultation. The MIMU (Opal sensors, APDM, Portland, OR, USA) was aligned with the set square by means of a custom adapter. Sensor data (triaxial accelerometers, gyroscopes, and magnetometers) were collected at 128 Hz.
Time synchronization between both systems was ensured by different procedures. During the experiment, an external trigger executed by the launch of the MIMU system measure triggered the optoelectronic system measure. However, after resampling the measurement from the optoelectronic system by cubic interpolation at 128 Hz to be consistent with the inertial measurement, postprocessing showed a time-shift between the data. Therefore, another procedure, depicted in Annex A, was applied to time-synchronize the data. The set square was moved manually in the calibrated volume of the optoelectronic system. In order to analyze the effect of external accelerations experienced by the MIMU on covariance matrices, three different intensity levels of movement were performed on the set square. To define these three intensity levels, the set square was moved as rapidly as possible and the corresponding values were recorded. The three intensity levels were then defined by creating three equal intervals between zero and the maximal recorded values. These levels were designated by the adjectives slow, intermediate, and fast. Table 1 lists the main kinematic characteristics corresponding to each movement intensity level. Given the likely presence of static or quasi-static periods during human motion, we subsequently added a static acquisition to our protocol to analyze this particular situation. Each movement was executed for a duration of 10 min so that the possible drift caused by the gyroscopes had time to intervene [18]. The movements were also initiated by a short static phase to estimate the initial bias of the gyroscopes, as well as the initial orientation of the set square. The nature of the imposed movements was, to the extent possible, "random" and of "globally constant" intensity. The global acceleration and angular velocity were controlled retrospectively. Five recordings for which the kinematics did not follow the same movement intensity level over the 10 min period were performed again. The movements were not collected during typical human movement, such as gait, because this could have resulted in reduced combinations of acceleration, angular velocity, and magnetic field values.

Reflective markers
Reference ⃗ ⃗ ⃗ Figure 1. Set square equipped with both an MIMU and five reflective markers. The reflective markers can be used to define local coordinate systems associated with the set square.
The set square was moved manually in the calibrated volume of the optoelectronic system. In order to analyze the effect of external accelerations experienced by the MIMU on covariance matrices, three different intensity levels of movement were performed on the set square. To define these three intensity levels, the set square was moved as rapidly as possible and the corresponding values were recorded. The three intensity levels were then defined by creating three equal intervals between zero and the maximal recorded values. These levels were designated by the adjectives slow, intermediate, and fast. Table 1 lists the main kinematic characteristics corresponding to each movement intensity level. Given the likely presence of static or quasi-static periods during human motion, we subsequently added a static acquisition to our protocol to analyze this particular situation. Each movement was executed for a duration of 10 min so that the possible drift caused by the gyroscopes had time to intervene [18]. The movements were also initiated by a short static phase to estimate the initial bias of the gyroscopes, as well as the initial orientation of the set square. The nature of the imposed movements was, to the extent possible, "random" and of "globally constant" intensity. The global acceleration and angular velocity were controlled retrospectively. Five recordings for which the kinematics did not follow the same movement intensity level over the 10 min period were performed again. The movements were not collected during typical human movement, such as gait, because this could have resulted in reduced combinations of acceleration, angular velocity, and magnetic field values. To verify the repeatability of the parameters of the covariance matrices identified, the three movements were repeated three times (on three separate days) for each of the three Opal APDM MIMUs tested. Therefore, a total of 27 motion capture sessions were made for these three sensors.

Gold Standard Measure
Orientation of the set square was estimated using the QUEST algorithm [19], which gives an optimal solution in a least-squares sense, to the Wahba problem [20]. The Wahba problem seeks to find a unique quaternion q that transforms all of the vectors that define the local coordinate system, based on the reflective markers seen in Figure 1, into the global coordinate system of the optoelectronic system. This algorithm makes it possible to consider all five markers for orientation, which tends to reduce measurement errors made on each of the markers. According to the QUEST algorithm, the Wahba's problem can be formulated as: where K is the matrix defined by: and B = ∑ k w k v T k w k and v k represent the vectors defined in the local and global coordinate systems respectively. Since three markers were placed on each part of the set square, three possibilities exist to define the x and the y axis of the local coordinate system (see Figure 1 for the representation of the x and y axis). Therefore, 2 × 9 = 18 local coordinate systems can be obtained, since it is possible to define the local coordinate system by keeping the x axis (such that . The same importance (and therefore the same weight in the QUEST algorithm) was given to each defined local coordinate system. The optimal quaternion corresponds to the associated normalized eigenvector, which pertains to the largest positive eigenvalues of K [19].
The orientation error resulting from the measured position of the reflective markers (mainly white noise) was then evaluated. First, the Allan Variance method [21] was applied to the five reflective markers during a static acquisition. For the five reflective markers, the white noise characteristic was of 5.7 µm in mean. The impact of this noise over the orientation was estimated by entering the identified white noise into simulated data that were used as input of the QUEST algorithm. This simulation enabled us to evaluate this gold-standard measurement error to 3 × 10 −30 on the set square orientation estimate.

MIMU Data Processing
The procedure was applied to three MIMUs (Opal sensors, APDM, Portland, OR, USA). Sensor data (triaxial accelerometers, gyroscopes, and magnetometers) were collected at 128 Hz and fused in a custom Kalman filter, as seen in Figure 2. According to the manufacturer, the noise density was 1.3 mm/s 2 / √ Hz for the accelerometer, 2.2 mrad/s/ √ Hz for the gyroscope, and 160 nT/ √ Hz for the magnetometer.
As the measurement model that relates sensor data to the orientation (expressed as a quaternion) is non-linear, we built an extended Kalman filter (EKF). To facilitate quaternion computation and to lighten the state vector, we used a multiplicative indirect Kalman structure [22]. The state vector was augmented with gyroscope bias, the behavior of which was modeled as a random walk. Angular velocity was thus corrected before the numerical integration. More details about this Kalman filter can be found in Annex B. lighten the state vector, we used a multiplicative indirect Kalman structure [22]. The state vector was augmented with gyroscope bias, the behavior of which was modeled as a random walk. Angular velocity was thus corrected before the numerical integration. More details about this Kalman filter can be found in Annex B. The reference orientation of the set square obtained by the optoelectronic system was expressed in the global coordinate system of the optoelectronic system, whereas the orientation resulting from the MIMU measurement was expressed in the north-east-down (NED) geographical coordinate system (the first axis defines the northern axis, the second the eastern axis, and the last, the descending vertical). For the purpose of comparison, the orientations measured by the two systems must therefore be expressed in the same coordinate system. To this end, we identified the quaternion that represents the transformation of the optoelectronic system coordinate system into the NED coordinate system.
To identify this transformation quaternion, we exploited the measurements obtained during the experimental procedure. The identification of this rotation was carried out in two steps. Firstly, the roll-pitch components of the transformation were identified at the beginning of the acquisition during a static phase, by taking the average quaternion [23] defined as follows: → * ⊗ Secondly, the yaw component of the transformation was identified, re-using the QUEST algorithm throughout the whole movement. This time, the optimal rotation being assessed corresponded to the transformation between the transformation → and the NED coordinate system. More details about the procedure is provided in Annex A.

Identification of the Kalman Parameters
As part of a multiplicative structure of the Kalman filter, the covariance matrix of the process noise was directly defined from the standard deviations characteristic of the noises of measurement [7]. The process and measurement and covariance matrices were then defined as follows: . . , . , and .  The reference orientation of the set square obtained by the optoelectronic system was expressed in the global coordinate system of the optoelectronic system, whereas the orientation resulting from the MIMU measurement was expressed in the north-east-down (NED) geographical coordinate system (the first axis defines the northern axis, the second the eastern axis, and the last, the descending vertical). For the purpose of comparison, the orientations measured by the two systems must therefore be expressed in the same coordinate system. To this end, we identified the quaternion that represents the transformation of the optoelectronic system coordinate system into the NED coordinate system.
To identify this transformation quaternion, we exploited the measurements obtained during the experimental procedure. The identification of this rotation was carried out in two steps. Firstly, the roll-pitch components of the transformation were identified at the beginning of the acquisition during a static phase, by taking the average quaternion [23] defined as follows: Secondly, the yaw component of the transformation was identified, re-using the QUEST algorithm throughout the whole movement. This time, the optimal rotation being assessed corresponded to the transformation between the transformation q roll−pitch opto→ned and the NED coordinate system. More details about the procedure is provided in Annex A.

Identification of the Kalman Parameters
As part of a multiplicative structure of the Kalman filter, the covariance matrix of the process noise was directly defined from the standard deviations characteristic of the noises of measurement [7]. The process Q and measurement R a and R m covariance matrices were then defined as follows: From this assumption, σ g , σ a , and σ m , denote respectively the standard deviation of the gyroscope, accelerometer, and magnetometer global errors (noises, external perturbations, calibration defects etc.), while σ bg represents the standard deviation of the gyroscope bias. The four parameters that were identified will be henceforth referred to as the Kalman parameters.
The process used to identify the four parameters is presented in Figure 3. For each set of Kalman parameters, the orientation was computed by means of the Kalman filter. On each occasion, this orientation was compared to the orientation resulting from the optoelectronic measurement. calibration defects etc.), while represents the standard deviation of the gyroscope bias. The four parameters that were identified will be henceforth referred to as the Kalman parameters.
The process used to identify the four parameters is presented in Figure 3. For each set of Kalman parameters, the orientation was computed by means of the Kalman filter. On each occasion, this orientation was compared to the orientation resulting from the optoelectronic measurement. This comparison enabled us to calculate the quaternion orientation error as follows:

⊗ *
This quaternion thus defined the rotation transformation from the orientation obtained with the optoelectronic system into the orientation obtained with the Kalman filter. Returning to the definition of the quaternion, we could then calculate the absolute angle corresponding to this rotation:

2.
Finally, the quality of orientation derived from the Kalman filter was characterized by the root mean squared error of this error angle over the 10 min of acquisition. This value, named RMSe, was obtained by the formula: 1 where N represents the number of measures over the 10 min.
To reduce the possible combinations of Kalman parameters, a set of ten values per parameter was chosen. This set of parameters was initially based on the characteristic values of measurement noises identified by the Allan's Variance method as seen in Table 2. These values were considered to be the minimal value for each Kalman parameter, because they should represent the minimum measurement errors. In fact, by definition they address only two stochastic disturbances (white noise and bias instability). So, to integrate the additional errors (other stochastic perturbations, calibration and model errors, and external perturbation), the nine other values of the set consisted of nine higher values. To explore a wide range of values, the ten values were regularly spaced on a logarithmic scale. This initial set of parameters was then refined to target the optimal parameters that were originally identified. We proposed to identify the parameters in two steps. First, the error was evaluated only in terms of roll-pitch, the effect of the magnetometers being therefore discarded. This approach made it possible to consider only three parameters initially: σg, σbg, and σa. At the end of this first step an This comparison enabled us to calculate the quaternion orientation error as follows: This quaternion thus defined the rotation transformation from the orientation obtained with the optoelectronic system into the orientation obtained with the Kalman filter. Returning to the definition of the quaternion, we could then calculate the absolute angle corresponding to this rotation: Finally, the quality of orientation derived from the Kalman filter was characterized by the root mean squared error of this error angle over the 10 min of acquisition. This value, named RMSe, was obtained by the formula: where N represents the number of measures over the 10 min.
To reduce the possible combinations of Kalman parameters, a set of ten values per parameter was chosen. This set of parameters was initially based on the characteristic values of measurement noises identified by the Allan's Variance method as seen in Table 2. These values were considered to be the minimal value for each Kalman parameter, because they should represent the minimum measurement errors. In fact, by definition they address only two stochastic disturbances (white noise and bias instability). So, to integrate the additional errors (other stochastic perturbations, calibration and model errors, and external perturbation), the nine other values of the set consisted of nine higher values. To explore a wide range of values, the ten values were regularly spaced on a logarithmic scale. This initial set of parameters was then refined to target the optimal parameters that were originally identified. We proposed to identify the parameters in two steps. First, the error was evaluated only in terms of roll-pitch, the effect of the magnetometers being therefore discarded. This approach made it possible to consider only three parameters initially: σ g , σ bg , and σ a . At the end of this first step an "ideal" value of σ bg was obtained. In the second step, the complete error was analyzed by integrating the effect of the magnetometers. We then considered the three parameters σ g , σ a , and σ m . Indeed, these three parameters are very closely linked and must be computed together.
This approach was adopted to provide a reasonable duration for the identification process. For example, given that the execution of the Kalman filter on a 10 min acquisition took about 30 s, the first identification step of three σ g , σ bg , and σ a lasted around 27 × 103 × 30 = 810,000 s, which already corresponded to nine days.

Optimality Criterion
For each movement, the average orientation errors on the three MIMUs tested, and on the three repetitions, were exploited. From the nine acquisitions at each intensity of movement, average and standard deviation of the RMSe were computed for each combination of Kalman parameters. To define the optimal values, it was proposed to consider a parameter defined by the sum of the RMSe and its associated standard deviation. Our optimality criterion was then to minimize this parameter. As such, this parameter takes into account not only the orientation error but also its variability during the identification process, which is more suitable for generalizing these results across more than three sensors.

Validation Procedure
To test the efficiency of the Kalman filter based on the optimal parameters identified, we performed a new experiment. First, we imposed on an MIMU from the same manufacturer (Opal sensors, APDM, Portland, OR, USA), that was not used during the identification procedure, a succession of movements, varying in intensity for a duration of 10 min. Then, we obtained orientations by setting the Kalman filter in five different ways. We first set the Kalman filter with the values of measurement noise identified by the Allan variance method, as seen in Table 2, then with the optimal Kalman parameters obtained for slow, intermediate, and fast movements. In addition, we implemented an adaptive Kalman filter as proposed by several authors [24,25]. In this study, we set the Kalman filter so that it selected automatically the set of parameters corresponding to the movement intensity, depending on the instantaneous acceleration and rotational speed. To prevent the Kalman filter from "oscillating" between two intensity modes, we defined a minimum duration of 10 s between each mode change. We will henceforth refer to this setting as "the adapted Kalman filter".
The error angle between the orientations obtained with each of the five different settings of the Kalman filter and the orientation obtained with the optoelectronic system was computed as well as the RMSe as previously depicted over the 10 min of acquisition.
The orientation error for the orientation obtained from the manufacturer algorithm was also computed. According to the manufacturer, their algorithm uses a state space model with a causal Kalman Filter. Different settings are possible. We chose the "Variable Weight Magnetometer" setting since according to the information provided by the manufacturer "it is the best choice for an environment similar to where the field calibration was last performed, and where having an absolute heading reference is desired", which was the case here.

Identification of the Kalman Parameters
Beyond identifying optimal values to assign to Kalman parameters σ g , σ bg , σ a , and σ m , the analysis of the results revealed interesting information concerning the behavior of the Kalman filter. For instance, Figure 4 represents, for the intermediate intensity movements, the RMSe translating orientation error between the orientation estimated with the optoelectronic system and that obtained with the Kalman filter, depending on the values of the Kalman parameters σ g and σ a . Most importantly, Figure 4 shows that a bad parameter choice can induce a huge error (up to 65 • for this movement), which confirms that appropriate selection of the Kalman parameters is essential.  The area on the extreme right of the figure corresponds to high values assigned to σg and low values assigned to σa. In other words, these values force the Kalman filter to give maximum importance to the measurement resulting from the accelerometers, and to neglect the measurement from the gyroscopes. Conversely, the area on the far left corresponds to low values assigned to σg and high values assigned to σa. In other words, these values force the Kalman filter to give maximum importance to measurement from gyroscopes, and to neglect the measurement resulting from the accelerometers. These two areas usually lead to significant errors. Indeed, the error observed in the right-hand area translates into an error in orientation estimation due to external accelerations experienced by the MIMUs, while the error observed in the left-hand area corresponds to the drift phenomenon resulting from the numerical integration of noisy angular velocity. The zone containing the optimal values is therefore located between these two particular areas, since it reflects a compromise between the gyroscope measurement and the measurement resulting from the accelerometers. Table 3 provides the optimal values for each Kalman parameter and the RMSe obtained for these values. These results confirm that the Kalman parameters are related to movement intensity. It is especially noted that σa is proportional to movement intensity (directly related with the external The area on the extreme right of the figure corresponds to high values assigned to σ g and low values assigned to σ a . In other words, these values force the Kalman filter to give maximum importance to the measurement resulting from the accelerometers, and to neglect the measurement from the gyroscopes. Conversely, the area on the far left corresponds to low values assigned to σ g and high values assigned to σ a . In other words, these values force the Kalman filter to give maximum importance to measurement from gyroscopes, and to neglect the measurement resulting from the accelerometers. These two areas usually lead to significant errors. Indeed, the error observed in the right-hand area translates into an error in orientation estimation due to external accelerations experienced by the MIMUs, while the error observed in the left-hand area corresponds to the drift phenomenon resulting from the numerical integration of noisy angular velocity. The zone containing the optimal values is therefore located between these two particular areas, since it reflects a compromise between the gyroscope measurement and the measurement resulting from the accelerometers. Table 3 provides the optimal values for each Kalman parameter and the RMSe obtained for these values. These results confirm that the Kalman parameters are related to movement intensity. It is especially noted that σ a is proportional to movement intensity (directly related with the external acceleration). These results also highlight the particular behavior of static condition because of which the Kalman filter must grant much less importance to gyroscopes (a hundred times less) in order to counteract bias instability. It can also be seen that the RMSe could not be reduced to under 13.7 • in fast movements. Figure 5 represents, over the whole movement, the orientation error computed between the orientation obtained with the optoelectronic system and those obtained with the four different settings of the Kalman filter based on the identified Kalman parameters. Because the orientation obtained with the Kalman filter based on the noise identified by the Allan's Variance method was particularly bad, for greater clarity we chose not to represent it in Figure 5. It can also be seen that the RMSe could not be reduced to under 13.7° in fast movements. Figure 5 represents, over the whole movement, the orientation error computed between the orientation obtained with the optoelectronic system and those obtained with the four different settings of the Kalman filter based on the identified Kalman parameters. Because the orientation obtained with the Kalman filter based on the noise identified by the Allan's Variance method was particularly bad, for greater clarity we chose not to represent it in Figure 5. The RMSe computed during the whole movement was 10° for the orientation obtained with the adapted Kalman parameters, and 46°, 16°, and 15° for the orientation obtained with the Kalman parameters identified respectively on slow, intermediate, and fast movements. The RMSe was 95° for that obtained with the Kalman parameters based on the noise identified by the Allan's Variance method. With the manufacturer algorithm, the RMSe computed during the whole movement was 20° for the orientation obtained, as seen in Figure 6. The RMSe computed during the whole movement was 10 • for the orientation obtained with the adapted Kalman parameters, and 46 • , 16 • , and 15 • for the orientation obtained with the Kalman parameters identified respectively on slow, intermediate, and fast movements. The RMSe was 95 • for that obtained with the Kalman parameters based on the noise identified by the Allan's Variance method. With the manufacturer algorithm, the RMSe computed during the whole movement was 20 • for the orientation obtained, as seen in Figure 6.

Results from the Validation Procedure
The RMSe computed during the whole movement was 10° for the orientation obtained with the adapted Kalman parameters, and 46°, 16°, and 15° for the orientation obtained with the Kalman parameters identified respectively on slow, intermediate, and fast movements. The RMSe was 95° for that obtained with the Kalman parameters based on the noise identified by the Allan's Variance method. With the manufacturer algorithm, the RMSe computed during the whole movement was 20° for the orientation obtained, as seen in Figure 6.

Discussion
In this study, a method for initializing covariance matrices has been proposed by comparing the inertial measurement with an optoelectronic measurement used as reference in the context of human motion analysis.
Results highlight the need to select the Kalman parameters required to define these covariance matrices rigorously, otherwise the performance of the Kalman filter can be seriously altered. This identification is even more important in that the optimal values obtained are relatively different from the quantities estimated by the Allan Variance method. Moreover, with the Kalman filter set to the quantities estimated by the Allan Variance method, considerable orientation error was found. The main advantage of the proposed approach is that it is global. Therefore, errors taken into account gather stochastic errors as well as behaviors not taken into account by the calibration models, or by the modeling carried out within the Kalman filter.
Our results clearly show that the optimal Kalman parameter values depend on movement intensity, which is of great importance for applications requiring human movement characterization. We proposed in the present paper to build a Kalman filter such that it automatically selected the set of Kalman parameters corresponding to the movement intensity, depending on the instantaneous acceleration and rotational speed. With such an adaptation, during a 10 min movement, orientation RMSe was reduced to 10 • compared to the 15 • obtained with the Kalman filter set to the optimal Kalman parameters obtained for fast movements. The remaining 10 • of orientation RMSe can probably be mainly imputed to persistent problems that affect orientation estimation by MIMU for fast motion. Indeed, during the identification procedure, the difference between orientations obtained with the optoelectronic system and with our Kalman filter-even optimized-remained above 10 • for fast movements. In that case, the optimal behavior of the Kalman filter leads us to place much more importance on gyroscopes than on accelerometers. Therefore, the drift phenomenon resulting from the numerical integration of noisy angular velocity cannot be avoided. This underlines the need for less noisy gyroscopes [26]. We have to admit that MIMUs nowadays include gyroscopes of higher quality than those present in the MIMUs used in this study. Noise density and bias stability of the gyroscopes are in fact greatly reduced in the new generation of MIMUs. For instance, the gyroscope noise density announced by the manufacturer was 0.05 • /s/sqrt(Hz) for the version of the sensors used in the present study, whereas it is of 0.025 • /s/ √ Hz for the latest version of these sensors [27]. Therefore, better results can be expected in the future. Moreover, our protocol designed a fast motion with 4 g acceleration and 700 • /s angular velocity during 10 min. Such intensity is hardly to be expected during 10 min of human movement.
To test the benefits of the proposed method on MIMUs of different quality, and from a different manufacturer, we applied this methodology to another MIMU (MicroStrain 3DM-GX4-25).
The corresponding results are presented in Appendix A. For this tactical-grade MIMU, the orientation error was also notably reduced (from 0.8 to 0.4 • ) despite the already small orientation error obtained with the manufacturer algorithm also based on Kalman filter, as seen in Table A1. Therefore, we are confident that the present methodology can be beneficial to all sorts of MIMUs in the context of human movement analysis.
Furthermore, let us notice that it is not straightforward to compare orientation errors obtained in the present study with orientation errors mentioned in previous studies. Not only were MIMUs used in the studies from different manufacturers and, as such, have not identical technical features, but also movement characteristics used to validate the methodology might have been different in terms of duration, linear acceleration, and angular velocity. Unfortunately, in previous studies movement characteristics are seldom fully documented. In the present paper, we obtain an orientation error of 1.5 • and 2.9 • for movements performed at respectively slow and intermediate intensity movement. Our results are close to Sabatini's results who presented an orientation error of 1.63 • with a new extended Kalman Filter [8]. However, as this previous study did not provide information regarding the movement intensity level, additional interpretation of the differences between the results from Sabatini's study and the present one would remain speculative.
Our present study focused on Kalman filter. Recent publications showed complementary filters as a promising alternative for fusing inertial data to obtain orientation [10,11]. For instance, a recent paper has tested an improved explicit complementary filter that uses only acceleration and gyroscope data as inputs [28]. The orientation error for the yaw angle remains close to 1 • but the nature of the movement tested is here again not known. Further investigations should then be conducted in order to assess thoroughly the advantages and disadvantages of complementary filters comparatively to Kalman filters.
Our proposed method has the disadvantage of being painstaking to achieve, not only in terms of experimentation but also in terms of computation, since many combinations of Kalman parameters have to be tested. To speed up the process, an optimization algorithm could be implemented to run the Kalman filter a small number of times. Other studies of minimal movement duration should be also performed to reduce the experimentation time. The proposed methodology required an optoelectronic system to provide the reference orientation. Robotic devices can probably be used to substitute this optoelectronic system and also the tester who is required to move the set square.
In the present paper, the identification procedure has been proposed from the study of three different movement intensities. The behavior of the Kalman filter in these situations should be considered in any study of human movement that is likely to lead to intermediate or mixed intensity movements. It is possible at first to simply exploit the Kalman parameters corresponding to the closest intensity tested, as proposed in the validation section of the present study. However, to go further, a continuous approach to Kalman parameter management could be considered. Indeed, the analysis of the roll-pitch component of the error showed that the optimal behavior of the Kalman filter was obtained by respecting a certain ratio between σ a /σ b function of movement intensity. A procedure that defines the evolution of this ratio according to the measured accelerations could then be proposed. We have in this study considered Kalman parameters common to all three Opal APDM MIMUs tested. However, each MIMU is slightly different, such that the assignment of specific parameters to each sensor could be advantageous. Such a choice would, however, require the reproduction of this identification process for all MIMUs.
One can also notice that the movements performed were not "natural" human movements, such as walking. Therefore, the real benefits of the presented method for human movement analysis must be further probed.
Finally, only the influence of movement intensity on the Kalman parameters has been tested. Additional, similar studies could also be conducted to specify the influence of magnetic perturbations on the optimal Kalman parameters.

Conclusions
The method presented provides a unique solution to the problem of identifying the optimal initial covariance matrices values for Kalman filtering. With these optimized initial covariance matrices, the orientation estimation was greatly improved. However, the methodology should be improved to reduce the duration of the whole process.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Synchronization of the MIMUs and Optoelectronic Systems
Appendix A. 1

. Temporal Synchronization
A process of synchronization of the two systems was initially planned by managing the triggering of the optoelectronic system from a rising edge executed by the MIMU measurement system. However, the effectiveness of this method was not satisfactory, because a significant time shift remained between the two measures.
We then implemented an optimization process to minimize the gap between the angular velocity measured by the two systems with a least squares approach. For the optoelectronic system, the angular velocity expressed in the local coordinate of the set square was obtained by discrete time-derivation of the quaternion angle as follows:

. Temporal Synchronization
A process of synchronization of the two systems was initially planned by managing the triggering of the optoelectronic system from a rising edge executed by the MIMU measurement system. However, the effectiveness of this method was not satisfactory, because a significant time shift remained between the two measures.
We then implemented an optimization process to minimize the gap between the angular velocity measured by the two systems with a least squares approach. For the optoelectronic system, the angular velocity expressed in the local coordinate of the set square was obtained by discrete timederivation of the quaternion angle as follows: 2 and then ∆ ∆ Figure A1. Comparison of the angular velocity obtained from the MIMU (in red) with that obtained from the optoelectronic system (in black) at the beginning of the 10 min acquisition (at the top) and at end (at the bottom).
Nevertheless, even after making this correction, a time lag appeared at the end of the acquisition. As illustrated by Figure A1, the two curves were well synchronized at the beginning of the acquisition but were no longer so at the end of the acquisition (i.e., after 10 min of measurement). This systematic observation revealed a problem of temporal proportionality between the two measurement systems. Figure A1. Comparison of the angular velocity obtained from the MIMU (in red) with that obtained from the optoelectronic system (in black) at the beginning of the 10 min acquisition (at the top) and at end (at the bottom).
Nevertheless, even after making this correction, a time lag appeared at the end of the acquisition. As illustrated by Figure A1, the two curves were well synchronized at the beginning of the acquisition but were no longer so at the end of the acquisition (i.e., after 10 min of measurement). This systematic observation revealed a problem of temporal proportionality between the two measurement systems.
A new optimization process was then implemented to identify the temporal homothety to apply to the optoelectronic data. Once applied, the synchronization of data throughout the acquisition was obtained. According to this optimization process, the temporal proportionality coefficient to be applied was systematically 0.003% or 30 ppm. This value shows that the progressive desynchronization observed was very weak. Indeed, this corresponds to an offset of 18 ms after an acquisition of 10 min, which is then equivalent to just over two points on a signal sampled at 128 Hz.
We attributed a higher confidence to the MIMU sensors, because their data are time-stamped, which is not the case for the data from the optoelectronic acquisition. Indeed, for the optoelectronic system, the time stamps are just deduced from the sampling frequency, without any control of the actual measurement times. This hypothesis was confirmed during the equivalent analysis conducted with another MIMU, the MicroStrain 3DMGX4-25, which had the same result. This result implies that characterization of the time of the optoelectronic measurement does not seem to be without fault.

Appendix A.2. Identification of the Quaternion Representing the Transformation of the Optoelectronic System Coordinate System into the MIMU Coordinate System
This procedure was carried out in two steps. Firstly, the pitch and roll components were identified and then the yaw component of the transformation between the optoelectronic coordinate system and the MIMU coordinate system.
For that purpose, during a static phase during which the set square was lying on the floor, at each time-step, the transformation quaternion between the two coordinate systems was identified according to the formula: q roll−pitch opto→ned = q * opto ⊗ q MI MU Then, the average quaternion was computed according to the Markley et al. method [23]. One could notice that this transformation quaternion had mainly a roll component of 180 • . In fact, the MIMU coordinate system is the north-east-down (NED) coordinate system. Therefore, its z axis is vertical, pointing downward, whereas the z axis of the optoelectronic coordinate system, which is perpendicular to the calibration wand placed on the floor during the calibration procedure, is almost vertical, pointing upward. An intermediate coordinate system equivalent to the NED coordinate system in terms of pitch and roll, but different in terms of yaw, was then obtained.
During the movement, the quaternions obtained with the optoelectronic system were then expressed in this temporary coordinate system as follows: The same procedure was not used during the second step. Indeed, the yaw component obtained from inertial measure is the component that suffers the most from inaccuracy. Moreover, the static phase being short, we preferred to estimate the yaw component during the whole movement. For these different reasons, we used a QUEST algorithm and only the magnetic field measure.
The QUEST algorithm seeks the optimal rotation in the least squares sense [29]. Here, we expressed two vectors in the NED coordinate systems: The QUEST method is then used to find the attitude matrix A, such that: Ar 1 = s 1 and Ar 2 = s 2 The quaternions obtained through the QUEST method are then averaged to obtain q y opto→ned . Quaternions measured from the optoelectronic system and expressed in the NED coordinate system can then be obtained as follows: with ω the mean angular velocity defined by: We consider an indirect approach; the prediction is then performed with the error state equation defined by the orientation error δθ and the bias error ∆b as: The dynamic model is here given by: with F c and G c , v the process noise formulated as: [û×] is the preproduct matricial product of a vector u.
The corresponding discrete state transition matrix Φ is given by: The noise covariance matrix is defined by: which becomes in a discrete case: with Q 11 = σ 2 g .T.I 3 + σ 2 b g .(I 3 .  The propagation of the state covariance matrix is then: The covariance matrix of the state vector is also updated by: P a,k = P k|k−1 − K a .H a .P k|k−1 The second step proposes a yaw correction based on the magnetometer measure m: The state vector is only defined by: z m = m the covariance matrix by: R m = σ 2 m .I 3 and the measurement matrix by: whereẑ m is the theoretical magnetic field computed from the orientation previously predicted: The residual is then: Suh proposes to apply this update only to the quaternion and not to the gyroscope bias by simplifying the state covariance matrix: P m,k = P a,k(1:3,1:3) 0 3,6 0 6,3 0 6,6 After that, the residual covariance matrix can be computed as: S m = H m .P m,k H T m + R m To limit the effect of the correction only to the yaw component, Suh [25] proposed to modify the computation of the gain as following:  σbg (rad/s) σg (rad/s) σa (m/s²) σm (µT) Allan's variance 1 × 10 −4 1 × 10 −3 2 × 10 −2 0.6 Identification 0 2 × 10 −4 1 4 As for the APDM MIMUs, Figure A2 depicts that the optimal behavior of the Kalman filter was obtained by respecting a certain ratio between σa/σb. The values of the Kalman parameters identified with the methodology presented in the paper are presented in Table A1. Once again, one can notice that the parameters identified are notably different from those estimated with the Allan variance method. Table A2 presents the orientation error in terms of root mean squared error between the orientation obtained with the optoelectronic system and the orientation obtained (a) with the Figure A2. RMSe between the orientation estimated with the optoelectronic system and the Kalman filter depending on the values of Kalman parameters σ g and σ a for the slow intensity movement.
The values of the Kalman parameters identified with the methodology presented in the paper are presented in Table A1. Once again, one can notice that the parameters identified are notably different from those estimated with the Allan variance method. Table A2 presents the orientation error in terms of root mean squared error between the orientation obtained with the optoelectronic system and the orientation obtained (a) with the manufacturer algorithm (b) with our Kalman filter set with the parameters identified with the methodology described in the paper. Table A2. Orientation error for slow intensity movement.

RMSe (degree)
Manufacturer algorithm 0.8 Kalman filter set with identified parameters 0.4