Error-state Kalman filter for lower-limb kinematic estimation: Evaluation on a 3-body model

Human lower-limb kinematic measurements are critical for many applications including gait analysis, enhancing athletic performance, reducing or monitoring injury risk, augmenting warfighter performance, and monitoring elderly fall risk, among others. We present a new method to estimate lower-limb kinematics using an error-state Kalman filter that utilizes an array of body-worn inertial measurement units (IMUs) and four kinematic constraints. We evaluate the method on a simplified 3-body model of the lower limbs (pelvis and two legs) during walking using data from simulation and experiment. Evaluation on this 3-body model permits direct evaluation of the ErKF method without several confounding error sources from human subjects (e.g., soft tissue artefacts and determination of anatomical frames). RMS differences for the three estimated hip joint angles all remain below 0.2 degrees compared to simulation and 1.4 degrees compared to experimental optical motion capture (MOCAP). RMS differences for stride length and step width remain within 1% and 4%, respectively compared to simulation and 7% and 5%, respectively compared to experiment (MOCAP). The results are particularly important because they foretell future success in advancing this approach to more complex models for human movement. In particular, our future work aims to extend this approach to a 7-body model of the human lower limbs composed of the pelvis, thighs, shanks, and feet.


Introduction
Human lower-limb kinematic measurements are critical for many applications including gait analysis, enhancing athletic performance, reducing or monitoring injury risk, augmenting warfighter performance, and monitoring elderly fall risk, among others [1][2][3][4][5]. Historically, most research studies are constrained to laboratory environments where camera-based motion capture systems (MOCAP) are used to measure body segment kinematics. A major disadvantage of laboratory-based studies is that experimental constraints (e.g., limited capture volume, artificial environment, and observers) may alter how subjects perform tasks, making it difficult to extrapolate results to unconstrained (real-life) environments [6,7]. Some research studies require continuous monitoring of human kinematics, rendering laboratory-based methods ineffective [8][9][10]. In addition, camera-based systems are relatively expensive, require long setup times, and require trained researchers [1,7]. Collectively, the above realities of laboratory-based experiments may significantly limit the findings and benefits of the research. Thus, there remains a strong motivation to advance the use of wearable sensors to measure human kinematics outside of the lab environment. The most utilized sensors for mobile kinematic measurements are inertial measurement units (IMUs) which contain three-axis accelerometers and angular rate gyroscopes (with some designs also including magnetometers, GPS, barometers, or other sensors). The resulting measurements can be integrated and/or differentiated to estimate the kinematics of the body segments to which they are attached [11]. However, because of measurement noise and finite sampling rates, kinematic variables estimated via numerical integration are subject to integration drift errors. Consequently, accurate estimates of kinematics from (noisy) IMU measurements must also correct for drift errors [10][11][12][13][14].
In the context of estimating lower-limb human kinematics with IMUs, approaches that leverage known kinematic conditions and relationships demonstrate success in correcting integration drift errors in certain applications [15]. One well-known example is the zero-velocity update (ZUPT) method for computing three-dimensional trajectories of the feet [13,16]. This method uses the fact that, during human walking, the foot (and attached IMU) must be nearly still (zero-velocity) at some point during each stance phase to correct for drift in the estimated foot velocity. Other methods successfully estimate joint angles for single joints. For example, IMU-based knee joint angle algorithms capitalize on the fact that the knee often acts as a hinge joint [17,18].
In addition to single-segment or single-joint methods highlighted above, recent work for multi-segment or multi-joint systems shows progress towards describing the human lower limbs. While proprietary products for such analyses exist (e.g., Xsens MVN Link, Noraxon Ultium Motion), they incorporate unspecified assumptions (hence, unspecified limitations) which is especially pertinent because accuracy of IMU-based methods are often task-specific [15]. Independent validation studies of these systems confirm such task dependence and also reveal that accuracy varies significantly between specific joint angles [19][20][21]. These proprietary products also have significant cost and rely on product-specific IMU hardware or even specialized wearable suits. Thus, a significant need exists for validated and well-documented methods to advance future research and applications.
Several methods exist for estimating the kinematics of the human lower limbs using a 7-body representations of the human lower limbs constituting the feet, shanks, thighs and hip. Ahmadi et al. [22] utilize a ZUPT method to estimate ankle position trajectories and combine those with individual segment orientation estimates to yield estimated lower-limb kinematics for straight walking on level ground and stairs. Optimization ensures the joint angles conform to assumed ranges of motion. Results are validated via comparison with MOCAP measurements for short trials (six passes through a MOCAP volume) that may not fully expose the accumulation of (long-term) drift error. The results demonstrate strong correlations (R >0.94) for joint angles, but only those restricted to the sagittal plane. Teufl et al. [23] employ an iterated extended Kalman filter to estimate lower-limb kinematics and with root-meansquare (RMS) joint angle differences (all three axes) below 6 degrees relative to MOCAP measures. Additionally, their method estimates RMS stride length and step width differences of 0.04 and 0.03 meters, respectively, compared to MOCAP [24]. However, their algorithm assumes level-ground (to correct vertical drift and to identify zero-velocity update times), which renders it unsuitable for quantifying gait on general (unconstrained) terrain as often encountered outdoors. Collectively, the limitations of the studies reviewed above point to the need for a general algorithm that accurately estimates lower-limb kinematics over long trials (i.e., greater than five minutes) and without assumptions of terrain morphology. This paper presents a novel error-state Kalman filter (ErKF) method for estimating lowerlimb kinematics using data from wearable IMUs. We use this method to estimate three degree-of-freedom (DOF) joint angles, stride length and step width. The method presented here extends Sola's formulation of the ErKF for a single IMU [12] to a multi-IMU formulation and also incorporates biomechanical measurement models to correct integration drift errors. In contrast to [22], the method is effective for long integration times (i.e., long trials). In contrast to [23], the method does not embed a level-ground assumption. As a first step towards developing a complete (7-body) model of the human lower limbs, we first consider herein an approximate 3-body mechanical model. Doing so enables careful formulation and study of all key modeling steps but within the context of a simpler model. Additionally, evaluation on this well-characterized mechanical model permits direct evaluation of the ErKF method without the confounding error sources associated with human subjects including uncertainties in joint center locations, joint axes, sensor-to-segment alignment parameters, increased joint complexity, and soft tissue artefacts. Thus, evaluation of this novel ErKF method on a well-characterized 3-body mechanical model (a "walker") is a critical step towards extension to a full (7-body) model of the human lower limbs. We demonstrate the success of the method by comparison to two reference data sets. In the first comparison, estimated kinematic variables are compared to ground truth obtained by simulation. In the second comparison, estimated kinematic variables are compared to those measured by MOCAP using an engineered 3-body walker.

3-body lower-limb model
As a step towards estimating the lower-limb kinematics of a human (i.e., a 7-body model), we employ a novel ErKF method on a simplified (3-body) model of the lower-limbs for walking (i.e., a "walker"). We utilize the ErKF to estimate the poses (positions and orientations) of each IMU attached to each segment of the model. This simplified model embeds the key challenge to accurate lower-limb kinematic estimation from body-worn IMU data; namely, utilizing well-conceived measurement models to correct integration drift errors. The simplified model consists of a pelvis and two legs attached to the pelvis by hinge joints. The leg lengths and pelvic width are comparable to human anthropometrics (0.92 m and 0.39 m respectively). This model is simulated in OpenSim [25] (Fig 1A) and also fabricated for experiments ( Fig 1B). The OpenSim model is modified from OpenSim's "Dynamic Walking Challenge" example [26].
The model includes three IMUs with one near the center of mass of the pelvis (i.e., a sacrum-mounted IMU) and one at the distal end of each leg (i.e., foot-mounted IMUs). The IMUs are attached via a weld joint in the OpenSim model and via athletic tape in the experiment.

Error-state Kalman filter method (ErKF)
While traditional extended Kalman filter equations are written with respect to the states of the system directly, the ErKF equations are written with respect to the errors in these states which then correct the estimated states. The ErKF demonstrates superior performance over a traditional extended Kalman filter for similar applications for aircrafts and robots due to key advantages including linearity, singularity avoidance, and simplicity [12,[27][28][29]. Recent work also shows great promise for using ErKF formulations to improve the accuracy of joint angle estimates for biomechanical applications [30].
In this formulation, each segment is treated as a free body (i.e., possesses six DOF) which yields superior estimates of joint angles versus a minimal DOF kinematic chain, as previously illustrated for the human arm [31]. This algorithm requires knowledge of the positions of the joint centers and the directions of the anatomical axes for each segment in the attached IMU's reference frame. We refer to these relationships as the sensor-to-segment alignment and assume they remain constant throughout a trial. Through the process model of the ErKF, IMU data from each body segment is integrated to estimate the time-dependent pose of the IMU (thus the time-dependent pose of the segment). Errors are corrected through known kinematic constraints (e.g., the joints between segments) and kinematic states (e.g., if a segment is momentarily at rest) that are incorporated in the measurement model. We describe the filter below.
States and error-states. This multi-IMU formulation for the ErKF draws from and extends Sola's formulation for a single IMU [12] which also provides detailed derivations of relations used below. To reduce the size of the state, we do not estimate sensor biases or the gravitational acceleration. Instead, we assume sensor biases and the gravitational acceleration are well characterized (and constant) and these are estimated using IMU data during an intentional still period at the start of a trial. Therefore, the state for the j th IMU, x j , is the (10x1) vector where p j is the (3x1) position vector of the accelerometer within the IMU in a world (i.e. lab- fixed) frame, v j is the (3x1) velocity vector of the same point, and q j is the (4x1) quaternion rotation vector (using Hamiltonian convention) that relates the IMU sense axis frame (hereafter called the IMU frame) to the world frame. In particular, we define q as the rotation quaternion that transforms a vector in the body-fixed frame (y b ) to its components in the world frame (y w ) per where � denotes quaternion multiplication and q � denotes the quaternion inverse. The error state for the j th IMU, δx j , is the (9x1) vector where δp j and δv j denote errors in the position and velocity, respectively, and δθ j is the (threecomponent) attitude error vector (assumed to be small) defined such that the quaternion error δq j obeys where b is the unit vector in the direction of δθ j (i.e., the axis of rotation) and k�k is the Euclidean vector magnitude. The full state x is the concatenation of the states of all n IMUs in the system, namely x ¼ . .
Similarly, the full error state δx is : ð6Þ The error-state covariance matrix (associated with the full error state) is denoted by P. Process model. The prediction step of the ErKF uses the process model for each IMU where b x j denotes the prediction of x j , the additional subscript k denotes the k th time-step, and u j denotes the IMU data (acceleration and angular velocity). The state of each IMU at time-step k+1 is then predicted from the state and IMU data at the previous time-step where Δt is the sampling period of the IMU, R is the rotation matrix corresponding to q, g is the gravitational acceleration vector (in the world frame), a j is the acceleration measured by the j th IMU, and ω j is the angular rate measured by the j th IMU. Note that the predicted state for each IMU (Eq 8) is independent of the other IMUs (each IMU is treated as an independent six DOF rigid body). During the prediction step, the covariance matrix P is also estimated as follows. The Jacobian of the process model for the j th IMU at time-step k with respect to its error state vector, F xj,k , follows from and S{w} applies the Rodrigues' rotation formula on the vector w Here, φ is the scalar magnitude of w and s is the unit vector in the direction of w. Due to the independence of the IMUs in the prediction step (Eq 8), the Jacobian of the full system process model relative to the error state at time-step k, F x,k , is where blkdiag denotes the block diagonal matrix composition. The process noise covariance for the j th IMU, Q j , is where s 2 a and s 2 o are the noise variances for the acceleration and angular rate signals, respectively, with the values being obtained from the manufacturer or through experiment. Note that this matrix is assumed constant for each IMU and thus we do not denote a time step. Again, due to the independence of the IMUs in the prediction step, these matrices can be concatenated to form the full system process noise covariance matrix, Q, per Thus, the prediction update for P is The process model equations above detail predictions for the state and the error-state covariance. Note that the estimated error-state mean is not calculated because it is always zero throughout the process model because the error-state mean is initialized to zero and is reset to zero following any measurement update.
Measurement model. In the absence of any measurement during this time step we use the predicted state and error-state covariance as the best estimates at time-step k+1, namely However, when measurements are available during this time step, we apply corrections using the measurement model to improve the estimates as follows. The measurement model takes the functional form where z is the observed measurement, h(x) is the expected measurement represented as a function of the state x, and c is Gaussian white noise with covariance C. Specific measurement models follow below. For each, we linearize the measurement equation by defining the Jacobian H evaluated at x Consistent with [12], we use the chain rule to decompose H as where H x depends on the measurement model and X δx depends only on the estimated orientation. Next, the Kalman gain, K, and error-state mean, c dx are updated per The error-state mean for each IMU, c dx j , updates its respective nominal state per After the nominal state mean is updated, the error-state mean is reset to zero and the errorstate covariance is updated per where G k is the Jacobian of the error-state reset operation with respect to the error state at time-step k, defined as where G j;k ¼ � � : ð27Þ We note that the addition of the G terms in Eq 25 in the present formulation differs from the covariance measurement update in a traditional extended Kalman filter formulation to account for the error-state mean reset (to zero) after each measurement update. The above process of prediction and measurement updates (when available) repeats each time step.
Next, we present four measurement models. The first two pertain to known kinematic states of the body segments (e.g., when the IMU is still) and the second two pertain to known kinematic constraints (e.g., constraints imposed by the two joints). Note that in the case of multiple measurements during a time step, a batch measurement update is used (i.e., all measurements are stacked and processed together).
Measurement model 1: ZUPT correction. We leverage the fact that a foot will be momentarily at rest sometime during the stance phase during gaits that do not induce significant slipping. Thus we employ, a zero-velocity update (ZUPT) correction for estimating foot trajectories which accurately describes gaits at normal walking speeds [13,16] through fast walking and running speeds [32]. Within our framework, the associated measurement equation becomes where h ZUPT (x) is the expected measurement for the ZUPT correction and v IMU is the (3x1) velocity vector for a foot-mounted IMU. This expected foot velocity is compared to the (virtual) observed measurement of the foot velocity when the IMU (foot) is (momentarily) still. Measurement model 2: Gravitational tilt correction. We also leverage the fact that when an IMU is still, the accelerometer in the IMU measures only gravitational acceleration and therefore functions as an inclinometer, thus enabling a gravitational (tilt) correction for IMU orientation. This correction yields the measurement model where h tilt (x) is the expected measurement for the gravitational correction and R is the rotation matrix for the still IMU. Note that in Eq 30, it is assumed that gravity acts opposite the third world-frame component (i.e., in the "-z" direction); however, this equation can be easily modified to accommodate other world-frame definitions. This expected measurement is compared to the observed measurement of tilt with a being the IMU acceleration. Note that we compare the direction of unit vectors (h tilt and z tilt ) rather than the full acceleration vector to mitigate the effects of discrepancies in magnitude caused by the IMU not being exactly still or the effects of sensor noise and bias. Measurement model 3: Joint center correction. Next, the joint center between two adjacent limbs must be approximately at the same position as deduced from the positions and orientations of those limbs [31]. For IMUs on adjacent limbs 1 and 2, the measurement equation becomes where h jc (x) is the expected measurement for the joint center correction, the subscript i = 1,2 denotes IMU i , r i denotes the known position of the joint center from IMU i (and resolved in the IMU frame), and R i denotes the rotation matrix for IMU i . The (virtual) observed measurement for the joint center correction is Measurement model 4: Joint axis correction. Similar to the joint center correction above, at times the joint axis must be the same as deduced from the orientations of the adjacent limbs (IMUs). An example of this correction arises when the knee is predominantly acting like a hinge [17,18] and the flexion/extension axes of the thigh and shank must be aligned in the world frame. This can be generalized for any pair of adjacent limbs 1 and 2 per where h JA (x) is the expected measurement for the joint axis correction and e i is the aligned joint axis (unit vector) deduced from IMU i in the frame of IMU i . The (virtual) observed measurement for the joint axis correction is Evaluation of ErKF method using two reference data sets We evaluate the performance of the ErKF method using two sets of reference data, namely: 1) simulated IMU data for the simulated walker with associated simulated ground truth results and 2) experimental IMU data for the physical walker with associated MOCAP results. These two data sets allow us to evaluate the performance of the method with increasing levels of model complexity and uncertainty (e.g., knowledge of sensor noise characteristics, sensor-tosegment alignment) and without several confounding error sources from human subjects that affect both the ErKF method and reference MOCAP estimates including uncertainty in sensor-to-segment alignment parameters (e.g., relationships between sensor and anatomical coordinate systems, joint center locations), increased joint complexity, and soft tissue artefacts. To this end, we compare estimated and reference hip joint angles, stride length, and step width for the walker simulation and experiment (i.e., reference data sets 1 and 2 above). The reported hip joint angles mirror the International Society of Biomechanics convention for human subjects [33] with body axes defined such that the x-axis points anteriorly, the y-axis superiorly, and the z-axis to the right (aligned with hinge joint axis) when the model is in a neutral upright pose (as in Fig 1B). The definitions of stride length and step width parameters are consistent with [24] with the following minor modifications: 1) identified footfall instances are used in place of initial contact times, and 2) the IMU position is used in place of the heel position for the 3-body model. S1 Appendix details the methods for detecting footfalls and still periods. Further details for each reference data set are provided next. Reference data set 1: ErKF method estimates for walker compared to simulation. We first evaluate the performance of the ErKF method via simulation because simulation enables assessment of the ErKF method independent of many confounding factors associated with experimental data. For example, in the simulation we specify the 1) sensor error (e.g., bias, noise) parameters, 2) sensor-to-segment alignment parameters, 3) measurement times (e.g. when a foot is stationary, when a joint acts as a hinge), and 4) ground truth data for comparison. We first compute generalized coordinate trajectories for a straight-line walk for the three body segments. The gait consists of 200 identical strides with a mean speed and stride length of 0.33 m/s and 0.73 m, respectively. Stance and swing angular trajectories were chosen to have a waveform similar to the simplest walking model [34]. The gait also contains (0.1 second) still periods following each ground contact, permitting clear identification of times of zero-velocity of the "feet". The OpenSim model is then driven with these computed trajectories and the BodyKinematics analysis tool in OpenSim computes the virtual IMU poses with respect to a fixed lab-frame. Simulated IMU data (accelerations and angular rates) that are free of noise and bias are calculated by differentiating these poses at a sampling rate of 512 Hz (sampling rate of the IMUs used for Reference Data Set 2). Finally, real (i.e., noisy) IMU data is simulated by adding prescribed zero-mean Gaussian noise to this data. The accelerometer noise value is taken from the specification sheet for the commercial IMU (Opal, APDM, ±16 g and ±200g accelerometers, ±2000 deg/s gyro) used in the experiments. The gyro noise is that from the same specification sheet plus additional noise (10 deg/hr drift) to account for both bias instability and angular random walk. The joint center measurement noise comes from [23]. All noise values are summarized in Table 1.
Poses of each segment throughout the walk are estimated by employing the ErKF method with the simulated IMU data. Because the joints are constrained to be pure flexion/extension, the joint axis and joint center measurement corrections are applied at each time step with constant joint centers and joint axes. The segment poses follow directly from the IMU poses because the sensor-to-segment alignment is constant throughout the trial. Additionally, for this simulation, the sensor-to-segment alignment is known exactly.
Zero-velocity measurements are applied at identified footfalls while gravitational correction measurements are applied at each still period while using the measurement noise reported in Table 1. The accuracy of ErKF estimated hip joint angles, stride lengths, and step widths are compared to prescribed values from the original gait trajectories.
Reference data set 2: ErKF method estimates for walker compared to MOCAP. Next, we evaluate the performance of the method on the walker during overground walking gait. A marker-based motion-capture (MOCAP) system (Vicon, 18 Vero V2.2 cameras) tracks positions of reflective markers on the model at 100 Hz. Seven reflective markers are attached to each segment (four to define the primary axes and three additional markers, see Fig 1B). Positional estimates of the markers are filtered with a 4 th order low-pass Butterworth filter at 20 Hz. Additionally, the attached IMUs yield sampled acceleration and angular rate data at 512 Hz.
IMU poses on the segment (for sensor-to-segment alignment) are defined by three reflective markers attached to each IMU (Fig 1B). Body frame axes for the three segments are determined as described in Reference data set 1. Joint center locations are estimated at the center of the T-joints at the hips using the known dimensions of the model. A single still frame at the beginning of the trial is used to determine the positions of joint centers, IMUs, and reflective markers in the body-fixed frames and are assumed constant. MOCAP estimates of segment orientations are determined as follows using a published optimization method [35]. For each segment and time step, we record all pairwise positions between the markers on the segment and compare them to the same from the still frame. A MATLAB implementation of the aforementioned optimization method [36] is used to estimate segment orientation. Segment orientation estimates that yield (mean residual) marker positional errors exceeding 0.01 meters are eliminated as they indicate misidentified markers or significant marker positional error. Short time gaps (<0.05 seconds) in the orientation estimates are filled using linear interpolation for the resulting Euler angles ("unwrapped" to account for discontinuities) following which a 4 th order low-pass Butterworth filter (20 Hz cut-off) is applied for data smoothing. Rotational alignment between segment axes and their associated IMU sense axes is computed using a singular value decomposition procedure [37] comparing body-fixed IMU (from gyro measurement) to segment (from MOCAP segment orientation estimate) angular velocity vectors. During the trial, a researcher manually operates (walks) the model back and forth through the 4.5-meter capture volume for ten minutes. The average sacrum velocity during straight walking is 0.44 m/s. Because the model is restricted to pure flexion/extension of the hip, joint center measurements are applied at every time step. Additionally, joint axis measurements are applied at every time step recognizing that the body-fixed z-axes (flexion/extension axes) of the pelvis and legs are aligned. As with the simulation, zero-velocity measurements are applied at all footfall instances and gravitational tilt measurements applied at all still periods. Measurement noise parameters used in the ErKF for data set 2 are the same as for data set 1 (Table 1).

Reference data set 1: ErKF method estimates for walker compared to simulation
We open by comparing estimated outputs from the walker model (using simulated IMU data) to ground truth from the OpenSim simulation. Comparisons are made for all full strides excluding the first (transition) stride for each leg (198 strides for the right leg, 199 strides for the left leg). Joint angle estimates from the model are compared to ground truth values at each sample throughout the trial. Mean, standard deviation (SD), and root-mean-square (RMS) differences for the three joint angles are summarized in Table 2. RMS differences are specifically included for comparisons to results from [23,24]. Additionally, range of motion (ROM) is estimated and compared for each stride. Summary statistics for ROM differences are also reported in Table 2 for flexion/extension but not for internal/external rotation or abduction/adduction as their true values are constantly zero for this model. Fig 2 shows the differences in the three joint angles as functions of time for the right hip joint over this exemplary long (7 minute) trial (results similar for left hip). Importantly, the results reveal no observable drift in the joint angle differences with time (slopes of linear fits of the joint angle differences versus time remain below 0.1 deg/hr across all joint angles). By contrast, without any filter corrections, the differences can grow to up to 10 degrees due to drift over this same time interval. Fig 3 illustrates the estimated flexion/extension angle compared to ground truth during the gait cycle where time is normalized by gait cycle time (which begins and ends with the instances of identified footfalls). Shown are the average (solid line) and one standard deviation (shaded region) across all strides.
We also report the accuracy of the estimated stride length and step width as summarized in Table 3. The mean differences are less than 1% of the average values for both stride length (0.73 m) and step width (0.39 m), while the RMS differences remain within 1% and 4% respectively.
Finally, we report the accuracy of the estimated foot IMU trajectories for each leg for the duration of the stride cycle. Fig 4 compares the forward, lateral and vertical coordinates of the right (Fig 4A) and left (Fig 4B) foot IMU to ground truth. Both mean (solid line) and one standard deviation (shaded region) are illustrated as functions of time (normalized by gait cycle). Note significantly enlarged scales for lateral and vertical displacements. Additionally, note that these trajectories show the relative displacement of the IMU center and not the position of the ground contact point itself (refer to Fig 1A); thus, negative vertical displacement (Fig 4) does not necessarily represent penetrations of the ground.

Reference data set 2: ErKF method estimates for walker compared to MOCAP
Next, we compare estimated outputs from the model using measured IMU data to those measured by MOCAP. The results below report the differences in estimated kinematical quantities

PLOS ONE
Error-state Kalman filter for lower-limb kinematic estimation: Evaluation on a 3-body model obtained by the two measurement modalities. Recall that the experimental procedure requires repeated walking through the MOCAP capture volume and thus consists of straight walks through the capture volume separated by sharp turns. Since the sharp turns do not represent human-like gait, we focus our evaluation only on the "straight walking strides". A straight walking stride is defined as one where the total displacement of the foot during the stride is no more than fifteen degrees from the primary direction of travel and the stride length is greater than 0.2 meters. We also exclude from analysis the transition stride following a turn for each leg.
Consider first the differences in joint angle estimates from the IMU-based and MOCAPbased methods using data from the straight walking strides (239 strides for the right leg, 215 strides for the left leg). Bland-Altman plots of Fig 5 [38] illustrate the 95% limits of agreement (mean difference ± 1.96 times the standard deviation) between the two measurement Table 3. Comparisons of stride length and step width estimates for simulation.
https://doi.org/10.1371/journal.pone.0249577.t003 modalities for the hip flexion/extension angle. Since the motion induces pure flexion/extension, we do not report Bland-Altman plots for internal/external rotation and abduction/adduction angles as the reference values are expected to be nominally zero. The mean differences for flexion/extension remain less than 1 degree for both hips with limits of agreement less than 3.2 degrees across both hips. The mean, SD, and RMS differences for all three joint angles and the ROM difference for flexion/extension are reported in Table 4. As with the Bland-Altman plots, we do not report range of motion results for internal/external rotation and abduction/adduction angles as the reference values are expected to be nominally zero.
Next, we evaluate how the differences in estimated joint angles vary with time over the entire ten-minute trial. Fig 6 illustrates the right hip joint angle differences versus time for all straight walking strides (similar results for left hip). While very small biases between the two joint angle estimates exist, the results reveal no observable drift in the differences over the tenminute trial (slopes of linear fits of the joint angle errors versus time remain below 1.8 deg/hr across all joint angles). By contrast, without any filter corrections, the differences can grow to up to 13 degrees due to drift over the 10-minute trial.
Consider next a comparison of the flexion/extension angle estimates through the gait cycle as reported in Fig 7, following the same procedure described above in the context of Fig 3. Illustrated are the mean (solid curves) and one standard deviation from the mean (shaded regions) for both measurement modalities. The largest differences in the means arise during the stance phase. However, the measured stride-to-stride variability in flexion/extension angle (the average width of the shaded regions) is only 0.2 degrees larger for the IMU-based estimates versus MOCAP estimates.  Finally, we assess the differences in estimated stride length, step width and the trajectories of both feet for all straight walking strides. The differences in estimated stride length and step width over all strides are reported in Table 5. The mean differences remain below 2% of the average value for both stride length (0.77 m) and step width (0.38 m), while the RMS differences remain below 7% and 5%, respectively. Fig 8 compares the forward, lateral and vertical coordinates of the right ( Fig 8A) and left ( Fig 8B) foot IMU to those measured by MOCAP. Both mean (solid line) and one standard deviation (shaded region) are illustrated as functions of time (normalized by gait cycle). Note significantly enlarged scales for lateral and vertical displacements. Additionally, note that these

PLOS ONE
Error-state Kalman filter for lower-limb kinematic estimation: Evaluation on a 3-body model trajectories show the relative displacement of the IMU center and not the position of the ground contact point itself (refer to Fig 1B); thus, negative vertical displacement (Fig 8) does not necessarily represent penetrations of the ground.

Discussion
This paper presents an IMU-based method that accurately estimates the kinematics of a simplified 3-body model of the human lower limbs for overground walking. The estimation method, developed using an error-state Kalman filter, fuses acceleration and angular rate data from three independent IMUs (one per rigid body) using four kinematic constraints. The kinematic constraints capture 1) foot zero-velocity updates, 2) gravitational tilt corrections, 3) joint center corrections and 4) joint axis corrections. The model is tested using two sets of comparison data, namely: 1) simulated IMU data from a simulated walker that yields ground truth results, and 2) experimental IMU data from a physical walker with associated MOCAP results. Results using simulated IMU comparison data demonstrate the success of the underlying IMU-based estimation model when all required inputs (e.g., noise parameters, sensor-to-segment alignments) are known and when critical measurement times (e.g., zero-velocity times) are also known. Under these conditions, the hip angle estimates remain highly accurate (RMS difference below 0.2 degrees) and without detectable drift error despite long-duration trials (~7 min). Additionally, stride length and step width estimates exhibit very small differences (RMS less than 1% and 4% respectively) when compared to ground truth values.
The experimental results for the second comparison data set with an experimental walker demonstrate excellent agreement between the IMU and MOCAP-based estimates for the hip joint angles, stride length, step width, and foot trajectories. In particular, the IMU and MOCAP-based estimates of the hip joint angles exhibit limits of agreement less than 3.2 degrees for the flexion/extension angle with RMS differences less than 1.4 degrees across all three hip angle axes. Additionally, RMS differences for stride length and step width remain below 7% and 5% respectively, compared to nominal values. Importantly, the differences in all kinematic variables did not appear to drift despite the longer duration trial (~10 minutes). Because the differences neither drift nor increase with time, there is every reason to anticipate similarly tight error bounds over even longer time periods. Moreover, the IMU-based estimates exhibit similar variation in the flexion/extension angle compared to MOCAP estimates; refer to the similar spread of the shaded regions in Fig 7. In the absence of a prior study like ours on a mechanical walker, we instead compare this study to prior work evaluating IMUbased methods on humans, while also recognizing the limitations of such comparisons. Recall that Teufl et al. [23] estimate the joint angles of the human lower limbs from IMU data (via an iterated extended Kalman filter) during a 6-minute walking test for comparison to MOCAP results [23]. Consequently, results from this study might be compared to those of [23] by focusing on the knee joint in [23] as it acts predominantly like a hinge during walking (as does the hip joint in the 3-body model herein). Doing so reveals that the RMS differences (relative to MOCAP) reported in [23] (~1.5 degrees) is remarkably similar to that reported herein (1.37 degrees), although also acknowledging that this comparison is limited due to obvious differences between the 3-body model for a walker and a 7-body model for a human. Importantly though, the model herein removes limitations in [23] including the assumption of level ground. Consequently, the method herein may hold great promise in extending to a 7-body model for a human and particularly for applications where the level ground assumption does not hold, including during activities of daily living and for outdoor sports and exercise. Accurate results for this second comparison data set may also be affected by uncertainties in model inputs such as sensor-to-segment alignment and footfall and still period detection. We mitigate these uncertainties by using MOCAP data to establish sensor-to-segment alignment and to manually correct misidentified footfalls. While this yields a method that is not truly "MOCAP-free", these topics (sensor-to-segment alignment and still period detection from IMU data alone) are themselves active areas of research for human applications [39][40][41].
The success of the novel ErKF method on a simplified mechanical model of the human lower limbs demonstrated in this study motivate its extension to a full (7-body) model. Extension to the full model for the human lower limbs presents many additional challenges including those due to complex human (versus mechanical) joints, uncertainties in model inputs (e.g., sensor-to-segment alignment, footfall detection), and soft tissue artefacts. Thus, future work must ensure that extensions of this ErKF method maintain accurate kinematic estimates despite these additional error sources. Importantly, future extensions should be evaluated against many different types of human gait, including abnormal gaits (e.g., due to injury or disease), to demonstrate its utility in a variety of clinical and other biomechanical applications. We note that the ErKF method developed here yields great promise for accurately estimating human lower-limb kinematics in such applications because it relies only on kinematic constraints that are largely independent of movement type.