Estimating Lower Limb Kinematics Using a Lie Group Constrained Extended Kalman Filter with a Reduced Wearable IMU Count and Distance Measurements †

Tracking the kinematics of human movement usually requires the use of equipment that constrains the user within a room (e.g., optical motion capture systems), or requires the use of a conspicuous body-worn measurement system (e.g., inertial measurement units (IMUs) attached to each body segment). This paper presents a novel Lie group constrained extended Kalman filter to estimate lower limb kinematics using IMU and inter-IMU distance measurements in a reduced sensor count configuration. The algorithm iterates through the prediction (kinematic equations), measurement (pelvis height assumption/inter-IMU distance measurements, zero velocity update for feet/ankles, flat-floor assumption for feet/ankles, and covariance limiter), and constraint update (formulation of hinged knee joints and ball-and-socket hip joints). The knee and hip joint angle root-mean-square errors in the sagittal plane for straight walking were 7.6±2.6∘ and 6.6±2.7∘, respectively, while the correlation coefficients were 0.95±0.03 and 0.87±0.16, respectively. Furthermore, experiments using simulated inter-IMU distance measurements show that performance improved substantially for dynamic movements, even at large noise levels (σ=0.2 m). However, further validation is recommended with actual distance measurement sensors, such as ultra-wideband ranging sensors.


Introduction
Human pose estimation is the tracking of position and orientation (i.e., pose) of body segments. Joint angles can then be calculated from the relative pose of linked body segments. Applications exist in robotics, virtual reality, animation, and healthcare (e.g., gait analysis). Traditionally, human pose is captured within a laboratory setting using optical motion capture (OMC) systems with up to millimeter level position accuracy [1] when properly configured and calibrated. However, the setup for OMC systems is time consuming and inconvenient (e.g., multiple markers are taped to the body) and requires considerable expertise. Recent miniaturization of inertial measurements units (IMUs) has paved the path toward inertial motion capture (IMC) systems suitable for prolonged use outside of the laboratory. Some examples of clinical applications in the current literature include determining level of spinal damage [2] and Parkinson's disease diagnosis [3]. Furthermore, developing a comfortable IMC for routine daily use may facilitate interactive rehabilitation [4,5], and allow the study of movement disorder progression to enable predictive diagnostics.
Commercial IMCs attach one sensor per body segment (OSPS) [6], which may be considered too cumbersome and expensive for routine daily use by a consumer due to the number of IMUs required. Similar works also exist in the literature using quaternions [7], Kalman filters (KF) [8], and particle filters [9]. Each IMU typically tracks the orientation of the attached body segment using an orientation estimation algorithm (e.g., [10,11]), which is then connected via linked kinematic chain, usually rooted at the pelvis. A reduced-sensor-count (RSC) configuration, where IMUs are placed on a subset of body segments, can improve user comfort, reduce setup time and system cost. However, using fewer wearable sensor units necessarily reduces the kinematic information available, which must otherwise be inferred from (i) our knowledge of human movement (e.g., enforcing mechanical joint constraints or making dynamic balance assumptions), or (ii) by using additional sensing modalities within each wearable sensor unit. Each approach will be described in the next subsections.

Existing Algorithms for Human Motion Tracking Using Fewer IMUs
The performance of algorithms with RSC configuration depends (i) on how the kinematic information of uninstrumented body segments is inferred and (ii) on how body pose is represented.
The kinematic information of body segments which do not have sensors attached to them may be inferred by the algorithm, either through data obtained in the past (i.e., observed correlations between co-movement of different body segments or data-driven approaches) or by using a simplified kinematic model of the human body (i.e., model-based approaches). Data-driven approaches (e.g., nearest-neighbor search [12] and bi-directional recurrent neural network [13]) are able to recreate realistic motion suitable for animation-related applications. However, these approaches are normally biased toward motions already contained in the database, which may limit their use in monitoring pathological gait. Model-based approaches reconstruct human movement using kinematic and biomechanical models (e.g., linear regression [14], constrained KF [15], extended KF (EKF) [16], particle filter [9], and window-based optimization [17]). The use of optimization-based estimators is sometimes favoured due to its relative ease to setup and ease of understanding. However, the algorithm can be very inefficient when tracking a larger number of dimensions (e.g., when tracking body pose over a long duration). To significantly increase the efficiency of the algorithm when estimating body pose across time, a recursive estimator can take advantage of the state substructure and reduce the state dimension being tracked [18].
Body poses are usually represented using Euler angles or quaternions [9,16]. However, recent work on pose estimation has shown that using a Lie group to represent the states of recursive estimator is a promising approach. Such algorithms typically represent the body pose as a chain of linked segments using matrix Lie groups to represent the orientation or pose of each body segment; specifically the special orthogonal group, SO(n), and special Euclidean group, SE(n), where n = 2, 3, are the spatial dimensions for humam body kinematics problems. The early works of Wang et al. [19] and Barfoot et al. [20] investigated representations and propagation of pose uncertainty, the former in the context of manipulator kinematics and the latter focused on SE(3), followed by the formulation of recursive estimators using Lie group representation (e.g., EKF [21] and unscented KF (UKF) [22]). Recent literature has reported the use of Lie group based recursive estimators to estimate human movement. Cesic et al. estimated the full body pose from OMC marker measurements and achieved significant improvements compared to an Euler angle representation [23]; and even supplemented the approach with an observability analysis [24]. Joukov et al. represented the human body using SO(n), tracking the pose using measurements from IMUs under an OSPS configuration [25]. Joukov et al.'s algorithm was tested using an arm tracking experiment, where the results improved especially during arm poses that cause a singularity when using an Euler angle representation (i.e., Lie group representation is singularity free).

Improving Human Pose Estimation Using Additional Sensor Measurements
Another approach is to supplement kinematic information from the IMU with another kind of sensor, which inherently increases cost and reduces battery life. Note that we will focus on systems that supplement pose estimation, not on the global position estimation of the subject (e.g., [26]). For example, IMCs can be supplemented with standard video cameras (e.g., fused using an optimization-based algorithm [27], and deep neural networks [28]) or depth cameras [29] at fixed locations in the capture environment, external to the subject. The combination of IMCs and portable cameras solves a weakness of OMCs (i.e., marker or body segment occlusion) and a weakness of IMCs (i.e., global position drift). However, the system still requires an external sensor that is carried by another person or requires some quick setup. IMCs can also be supplemented by distance measurements (using ultrasonic devices and KF in OSPS configuration [30], using constrained KF in RSC configuration [31]), removing dependence on any external non-body-worn sensor.

Novelty
This paper describes a novel human pose estimator that represents the state using Lie groups with the state propagated iteratively using a constrained EKF (CEKF) to estimate lower body kinematics for an RSC configuration of IMUs and inter-IMU distance measurements; the Lie group framework and inclusion of inter-IMU distance measurements, along with the exploration of its effect on pose estimation accuracy, are the major advancements made in this paper. It extends the work of [32] and builds on prior work of [15,31], but instead represents the state variables as elements of Lie groups, specifically SE(3), to track both position and orientation (whereas [15] only tracks position and assumes orientation measurements are noise-free). Furthermore, this paper presents in detail a novel Lie group formulation for assumptions typically used in human pose estimation (e.g., zero velocity update, constant body segment lengths, and a hinged knee joint). While not our focus here, our algorithm, with its SE(3) representation, is able to track the global position of body segments, taking into account IMU measurements during the prediction step, and a simpler implementation of certain measurement assumptions (e.g., zero velocity update), though at the expense of having an additional constraint step to ensure satisfaction of biomechanical constraints. The algorithm design was motivated by the need for a body pose representation that more closely models the human biomechanical system (without a dramatic increase in the dimensions of the tracked state) from which the missing kinematic information of uninstrumented body segments are inferred. The contributions of this paper advance the development of gait assessment tools for comfortable and long-term monitoring of lower body movement.

Mathematical Background: Lie Group and Lie Algebra
The matrix Lie group G is a group of n × n matrices (e.g., SE (3)). Mathematically speaking, it is also a smooth manifold with smooth group composition and inversion (i.e., matrix multiplication and inversion). The Lie algebra g can be represented in the vector space and is closely related to Lie group G. It represents a tangent space of a group at the identity element [33]. The elegance of Lie theory lies in its ability to represent pose using a vector space (e.g., Lie group G is represented by g) without additional constraints (e.g., without requiring R T R = I which is using a rotation matrix representation of orientation, or || q|| = 1 which is using a quaternion representation of orientation) [34].
The matrix exponential exp G : g→G (Equation (1)) and matrix logarithm log G : G→g relates (i.e., local diffeomorphism) the Lie group G and its Lie algebra g. The Lie algebra g is a n × n matrix that can be represented compactly in an n-dimensional vector space. A linear isomorphism (i.e., one-to-one mapping) between g and R n is given by operators [ ] ∨ G : g→R n and [ ] ∧ G : R n →g, which map between the compact and matrix representation of the Lie algebra g. Figure 1 shows an illustration of the said mappings.
Mapping between Lie group G, Lie algebra g, and a n-dimensional vector space.
[ Furthermore, the adjoint operator of a Lie group, Ad G (X), the adjoint operator of a Lie algebra, ad G (v), and the right jacobian, Φ G (v) (Equation (2)), where X ∈ G and [v] ∧ G ∈ g will be used in later sections. Multiplying an n-dimensional vector representation of a Lie algebra with Ad G (X) ∈ R n×n (i.e., the product Ad G (X)v) transforms the vector from one coordinate frame to another, similar to how rotation matrices transform points from one frame to another. ad G (v) is the Lie algebra of Ad G (X).A summary of the operators for Lie groups SO(3), SE(3), and R n will be explained in the next subsections. They will serve as building blocks to implement the algorithm being described by this paper. For a more detailed introduction to Lie groups refer to [18,34,35].

Special Orthogonal Group SO(3)
The special orthogonal group, SO(3) := R ∈ R 3×3 |RR T = 1, det R = 1 , represents orientation, where R is the typical 3 × 3 rotation matrix whose column vectors represent the x, y, and z basis vectors. The operations for SO(3) are listed below, and will serve as building blocks for SE(3), which will be described in the next subsection. Note that [a] ∧ SO(3) b is equivalent to the cross product of a and b. See Chapter 7 of [18] for details.
If φ/|φ| represents a unit vector axis we wish to rotate around, and |φ| is the angle by which we wish to rotate, then the rotation matrix, R, which will implement this rotation is given by Equation (4), which is also known as the Rodrigues' axis-angle rotation formula. When φ is very small, .
Furthermore, the Lie algebra adjoint, Lie group adjoint, and inverse operators are listed in Equation (5).

Special Euclidean Group, SE(3)
The special Euclidean group, SE(3) , represents orientation and translation, where T is the typical 4 × 4 transformation matrix, R is the rotation matrix, and t represents a coordinate point in Euclidean space. The operations for SE(3) are listed below. I i×i and 0 i×j denote i × i identity and i × j zero matrices. See Chapter 7 of [18] for details.
Lastly, we note the useful identity defined in Equation (11) (72)), which will be used to compute the Jacobians of our model later.

Real Numbers R n
In order to represent vector state variables (e.g., translation, velocity, and acceleration) and be consistent with how we used SE(3) to represent pose, we represented the real numbers s ∈ R n as SE(n) poses with position and no rotation. The operations for R n are listed below.
ad R n (s) = 0 n×n , Ad R n (S) = I n×n , Φ R n (s) = 0 n×n (14) Note that the multiplication of two elements of the Lie group (i.e., the exponential of s 1 and s 2 ) is equivalent to the vector addition of two elements of the Lie algebra (i.e., s 1 + s 2 ).

Algorithm Description
The proposed algorithm, L5S-3IMU, uses a similar model and assumptions to our prior works in [15,31], denoted as CKF-3IMU and CKF-3IMU+D, albeit expressed in Lie group representation. The algorithm uses measurements from three IMUs attached at the pelvis (sacrum) and shanks (slightly above the ankles), and the inter-IMU distance measurements to estimate the orientation of five body segments (i.e., pelvis, thighs, and shanks) with respect the world frame, W ( Figure 2). The Lie group representation enables the tracking of both position and orientation (note that CKF-3IMU only tracked position and assumed orientation measurements were noise-free). Figure 3 shows an overview of the proposed algorithm. L5S-3IMU predicts the ankle and pelvis positions through double integration of their linear 3D acceleration, and predicts the shank and pelvis orientation through integration of their linear 3D angular velocity. The orientation estimates are further updated using a third-party orientation estimation algorithm. Drift in the position estimates due to sensor noise, which accumulates quadratically in the double integration of acceleration, was mitigated through the following assumptions: (1) the pelvis position is approximated as the height of the pelvis when the leg(s) are unbent, or as informed by inter-IMU distance measurements, when available; and (2) the ankle velocity and height above the floor are zeroed whenever a footstep is detected. Furthermore, a pseudo-measurement equal to the current global position estimate of the pelvis and ankles is made with a fixed covariance to contain the ever-growing error covariance of the said states. Lastly, constant body segment lengths and hinged knee joints (one degree of freedom (DOF)) with limited range of motion (ROM) were enforced as biomechanical constraints. The pre-and post-processing parts remain exactly the same as the CKF-3IMU algorithm (e.g., acceleration due to body movement was calculated by expressing the acceleration of the instrumented body segment in the world frame using the orientation estimate and then subtracting acceleration due to gravity (i.e., g = [0 0 9.81] T ) [6]   Pre-processing calculates the body segment orientation, inertial body acceleration (calculated by expressing IMU acceleration with respect world frame using its current orientation estimate then subtracting gravity), and step detection from raw acceleration,ȃ k , angular velocity,ω k , and magnetic north heading,η k , measured by the IMU. The CEKF-based state estimation consists of a prediction (kinematic equation), measurement (orientation, pelvis height/inter-IMU distance measurement, covariance limiter, intermittent zero-velocity update, and flat-floor assumption), and constraint update (thigh length, hinge knee joint, and knee range of motion). Post-processing calculates the left and right thigh orientations, R R R lt k and R R R rt k .

System, Measurement, and Constraint Models
The system, measurement, and constraint models are presented below where k is the time step. X k ∈ G is the system state, an element of state Lie group G. Ω X k−1 , u k : G→R p is a non-linear function which describes how the model acts on the state and input, u k−1 , where p is the number of dimensions of the compact vector representation for Lie algebra g. n k is a zero-mean process noise vector with covariance matrix Q (i.e., where q is the number of dimensions of available measurements). D k ∈ G c is the constraint state, an element of constraint Lie group G c . c (X k ) : G→G c is the equality constraint function the state X k must satisfy (i.e., c (X k ) = D k ). Similar to [23,36], the state distribution of X k is assumed to be a concentrated Gaussian distribution on Lie groups (i.e., where µ k is the mean of X k and Lie algebra error ∼ N R p (0 p×1 , P)) [19].
The Lie group state variables X k model the position, orientation, and velocity of the three instrumented body segments (i.e., pelvis and shanks) as represents the pose (i.e., orientation and position) of body segment b relative to world frame W, and A v b is the velocity of body segment b relative to frame A. If frame A is not specified, assume reference to the world frame, W. The Lie algebra error is denoted as = [( T where the first three variables correspond to the Lie group in SE(3) while the latter three are for Refer to Sections 2.2 and 2.3 for definition of SE(3) and R n operators.

Lie Group Constrained EKF (LG-CEKF)
The a priori, a posteriori, and constrained state estimate (i.e., estimated mean of X k ) for time step k are denoted byμ − k ,μ + k , andμ + k , respectively. Note that the true state X k can be expressed as where µ k is one of the state means just mentioned with error, [ ] ∧ G . The a priori and a posteriori error covariance matrices are denoted as P − k and P + k , respectively. Note the error covariance is not updated at the constrain update step. The KF is based on the Lie group EKF, as defined in [36], where the state means (μ − k ,μ + k , andμ + k ) and state error covariance matrices (P − k and P + k ) are propagated by the KF at each time step.

Prediction Step
Prediction step estimates the a priori stateμ − k at the next time step. The mean propagation of the three instrumented body segments is governed by Equation (18) where Ω(μ + k−1 , u k ) (Equation (19)) is the motion model for the three instrumented body segments. Note that the state propagation may not necessarily respect the biomechanical constraints, so joints may become dislocated after this step. The input u k is defined in Equation (20), where the orientation and acceleration as obtained by the IMU attached to segment b with respect world frame W are denoted asȒ b k andȃ b k for b ∈ {p, ls, rs}, while the angular velocity as obtained by the IMU attached to segment The state error covariance matrix propagation is governed by Equation (21), where F k represents the matrix Lie group equivalent to the Jacobian of f (X k−1 , n k−1 ), Q is the covariance matrix of the process noise, and represents the linearization of the motion model with an infinitesimal perturbation . The process noise covariance matrix, Q, is calculated from the input-to-state matrix G and the noise variances of the measured acceleration and angular velocity, a and σ 2 ω , respectively.

Measurement Update
Measurement update estimates the state at the next time step by: (i) updating the orientation state using new orientation measurements of body segments from IMUs; by (ii) encouraging pelvis position to be above the feet, as informed by either some pseudo-measurement or inter-IMU distance measurements; and by (iii) enforcing ankle velocity to reach zero, and the ankle z position to be near the floor level, z f when step is detected. When only IMU measurements are available, (iia) pelvis z position is encouraged to be close to initial standing height, z p . When inter-IMU distance measurements are available, (iia) is not used. Instead, (iib) ankle distance is directly incorporated while pelvis position is inferred from inter-IMU distance measurements assuming hinged knee joints and constant body segment lengths. The a posteriori state meanμ + k is calculated following the Lie EKF equations below. (17) assuming m k = 0 and X k =μ − k , i.e., Z k = h(μ − k )) is the innovation/residual in Lie group G m brought to the vector representation of the Lie algebra space using the inverse exponential (i.e., logarithm) mapping.
H k can be seen as the matrix Lie group equivalent to the Jacobian of h(X k ), and is defined as the concatenation of H ori and H mp,k when inter-IMU distance measurement is not available. When inter-IMU distance measurement is available, H mp,k is replaced by H dist,k = [H T pla,k H T pra,k H T lra,k ] T . H ls,k and/or H rs,k are also concatenated to H k when the left and/or right foot contact (FC) is detected (See Equation (9) of [15]). Each component matrix will be described later. The measurement matrix Z k ∈ G m , measurement function h(X k ) ∈ G m , and measurement covariance noise R k are constructed similarly to H k , but combined using diag instead of concatenation (e.g., R k = diag(σ 2 ori , σ 2 mp )).
Orientation Update The orientation update utilizes the orientation measurement to update the state estimate as defined by Equation (30), with measurement noise variance σ 2 ori (9 × 1 vector).
H ori along with other components of H k are calculated by applying Equation (29) to their corresponding measurement functions, followed by tedious algebraic manipulation and first order linearization (i.e., exp([ ] ∧ ) ≈ I + [ ] ∧ ). The derivation for H ori (Equation (31)) can be solved trivially Pelvis Height Assumption The pelvis height assumption softly constrains the pelvis z position to be close to initial standing height z p as defined by Equation (32) (represented in vector space of its Lie algebra) and Equation (33), with measurement noise variance σ 2 mp (1 × 1 vector). This assumption is used only when inter-IMU distance measurement is not available. i x , i y , i z , and i 0 denote 4 × 1 vectors whose 1st to 4th rows, respectively, are 1, while the rest are 0; they are used below to select rows, columns, or elements from matrices. (36). Taking best estimate X k =μ − k gives us Equation (34).
Remember p T is a subvector of as defined in Section 3.1 and is the Lie algebra error of the state in its compact vector representation. Note that if measurement function h a (X k ) ∈ Lie group (15) and (13) (inverse of Lie group R n ). Finally, H mp,k is calculated as shown in Equation (36).
Zero Velocity Update and Flat Floor Assumption When step is detected, the ankle velocity is enforced to be zero and the ankle z position is brought to near the floor level, z f (i.e., flat floor assumptions). The corresponding measurement function is defined by Equation (37), with measurement noise variance σ 2 ls (4 × 1 vector). [ The zero velocity part of H ls,k (Equation (38)) and H rs,k can also be calculated trivially, while the flat floor assumption can be calculated similarly as H mp,k but the z position set to floor height, z f , instead of the pelvis standing height, z p .

Left and Right Ankle Distance Measurement
When the inter-IMU distance between the ankles,d lra k , is available, ankle distance measurement is incorporated as a soft distance constraint. The measurement function is defined by Equation (40), with measurement noise variance σ 2 lra (1 × 1 vector). τ lra (X k ) (Equation (39)) is the vector that points from the right ankle to the left ankle, where ls p la is the position of the left ankle expressed in left shank frame, and rs p ra is the position of the right ankle expressed in right shank frame. We have chosen that the ankles are at the origin of their respective shank frames. Note that matrix E converts homogeneous 4 × 1 coordinates to standard 3 × 1 coordinates (i.e., drops the 1 from the end of the 4 × 1 vector).
Pelvis-to-Ankle Distance Measurement In addition to the soft ankle distance constraint, the ankle to pelvis vector is inferred from the ankle to pelvis distance measurements while assuming hinged knee joints and constant body segment lengths. The measurement function is defined by Equation (45) The measurement pelvis to left ankle vector can be calculated from the measured pelvis to left ankle distance,d pla k as shown in Equation (46) which is the Lie Group reformulation of [31] (Equation (4)). In essence, Equation (47) calculates the most probably knee angle assuming hinged knee joint and constant body segment lengths, then Equation (46) adds the thigh (expressed in shank coordinate system with knee angleθ lk k ) and shank long axis to the hips to obtain the pelvis-to-ankle vector. See Appendix A for derivation. There are two solutions forθ lk k due to the inverse cosine in Equation (47). We chose theθ lk k value as that closer to the current left knee angle estimate from the prediction step. Note that this measurement function could also be formulated as a linearized Euclidean distance between the pelvis and ankle (i.e., similar to Equation (44)); however, a preliminary exploration of this approach showed poorer performance.
H pla,k is then calculated as shown in Equation (49). The right side of the pelvis-to-ankle distance measurement (i.e., h pra (μ − k ), Z pra , H pra,k ) can be solved similarly to the left side.
Covariance Limiter Lastly, the error covariance of the position estimates of the three instrumented body segments must be prevented from growing unbounded and/or becoming badly conditioned, as will occur naturally when tracking global position of objects without any global position reference. At this step, a pseudo-measurement equal to the current stateμ + k is used (implemented by Equation (50)) with some measurement noise of variance σ lim (9 × 1 vector). The covariance P + k is then calculated through Equations (51)-(53).

Satisfying Biomechanical Constraints
After the preceding updates, the joint positions or angles may be beyond their allowed range (i.e., knee hyperflexion). The constraint update corrects the kinematic state estimates to satisfy the biomechanical constraints of the human body by projecting the current a posteriori state estimatê µ + k onto the constraint surface, guided by our uncertainty in each state variable, which is encoded by P + k . The following biomechanical constraint equations are enforced: (i) estimated thigh long axis vector lengths equal the thigh lengths; (ii) both knees act as hinge joints (formulation similar to Section 2.3, Equation (4) of [9]); and (iii) the knee joint angle is within realistic range. The constraint functions are similar to Section II-E.3 of [15] but expressed under SE(3) state variables. The constrained stateμ + k can be calculated using the equations below, similar to the measurement update of [36] with zero noise, where C k = [C T L,k C T R,k ] T . C L,k is the concatenation of C ltl,k , C lkh,k , and C lkr,k ; the last matrix is not concatenated when the knee angle, α lk , is within its allowed range (i.e., α lk,min ≤ α lk ≤ α lk,max ). C ltl,k , C lkh,k , and C lkr,k corresponds to the biomechanical constraint for the left thigh length (ltl), left knee hinged joint (lkh), and left knee angle ROM (lkr), respectively, which will be described more later. C R,k can be derived similarly, while D k and c(μ + k )) are constructed similarly to Z k .
Thigh Length Constraint Firstly, the thigh length constraint is shown in Equation (59), where τ lt z (X k ) (Equation (58)) denotes the thigh long axis vector and d lt denotes the measured thigh length during calibration. p p lh is the position of the left hip expressed in pelvis frame, and ls p lk is the position of the left knee expressed in left shank frame. We have chosen that the left hip to be d p 2 to the left of the mid-pelvis origin, and the left knee to be d ls from the left shank origin (i.e., from the left ankle).
C ltl,k is calculated using Equation (60).
Following similar procedure to H lra,k , we obtain τ lt z (μ + k exp([ ] ∧ G )) = τ lt z (μ + k ) + Γ ltz (similar to Equation (42) Hinge Knee Joint Constraint Secondly, the hinge knee joint constraint as defined by Equation (62) is enforced by having the long (z) axis of the thigh to be perpendicular to the mediolateral axis (y) of the shank. For example, on the left leg, we would want r ls y be perpendicular to the thigh long axis vector τ lt z (μ + k ) (i.e., the dot product of r ls y and τ lt z (μ + k ) should be 0). Refer to Figure 2 for visualization. This formulation is similar to Section 2.3, Equation (4) of [9].

Knee Range of Motion Constraint
Thirdly, the knee ROM constraint is defined by Equation (69) and is only enforced if the knee angle, α lk , is outside the allowed ROM. The bounded knee angle, α lk , is calculated by Equation (67).

Equation (69) is obtained by expanding Equation (67) to Equation (68) which when rearranged gives us
[log(c lkr (X k ))] ∨ (i.e., Lie group representation of Equation (26) in [15]). Note that ls r lt z is the normalized thigh long axis expressed in the left shank frame. [log(c lkr (X k ))] ∨ = (E T ls ls r lt z =long axis of left thigh in shank frame Following a similar procedure to C lkh,k (i.e., replace i y in Equation (64) with ls r lt z ) and taking X k =μ + k , C lkr,k can be calculated from (70).

Post-Processing
The orientation of the pelvis and shanks are obtained from the stateμ + k .
The orientation of the right thigh, R rt+ , is calculated similarly.

Experiment
An extension of the dataset from [15] was used to evaluate our L5S based algorithms. The movements involved are listed in Table 1 (note the addition of dynamic movements), and were collected from from nine healthy subjects (7 men and 2 women, weight 63.0 ± 6.8 kg, height 1.70 ± 0.06 m, age 24.6 ± 3.9 years old), with no known gait abnormalities. Raw data were captured using a commercial IMC (i.e., Xsens Awinda) at 100 Hz sampling rate with IMUs attached to the pelvis and ankles, compared against a benchmark OMC (i.e., the setup followed Vicon Plug-in Gait protocol in a ∼4 × 4 m 2 capture area). The experiment was approved by the Human Research Ethics Board of the University of New South Wales (UNSW) with approval number HC180413. Frame alignment and yaw offset calibrations are similar to Section III-B of [15]. The experiments and algorithm were implemented using Matlab 2020a, with initial stateμ + 0 (i.e., position, orientation, and velocity) obtained from the OMC system (i.e., Vicon) and initial error covariance P + 0 set to 0.5I 27×27 . The variance parameters for the process and measurement error covariance matrix Q and R are shown in Table 2.  The inter-IMU distance measurements,d pla ,d pra , andd lra , were simulated by calculating the distance from the mid-pelvis to the left and right ankles and adding normally distributed positional noise with different standard deviations (i.e., σ dist ∈ {0, 0.01, . . . , 0.1, 0.15, 0.2} m). Each trial was simulated five times.
Lastly, the evaluation was done using the following metrics: (1) Mean position and orientation root-mean-square error (RMSE) (e.g., similar to [15,17] as shown in Equations (71) and (72)), where p b k and R b k are obtained from the benchmark OMC system,p b+ k andR b+ k are obtained from the algorithm. Note that as the global position of the estimate is still prone to drift due to the absence of an external global position reference, the root position of our system was set equal to that of the benchmark system (i.e., the mid-pelvis is placed at the origin in the world frame for all RMSE calculations). (2) Joint angles RMSE with bias removed (i.e., the mean difference between the angles over each entire trial was subtracted) and correlation coefficient (CC) of the hip in the sagittal (Y), frontal (X), and transverse (Z) planes and of the knee in the sagittal (Y) plane. Note that these joint angles are commonly used parameters in gait analysis. (3) Spatiotemporal gait parameters (e.g., total travelled distance (TTD) deviation, average stride length, and gait speed of the foot). Refer to Section III of [15] for more details. e pos,k = 1 N pos ∑ b∈DP || p b k −p b+ k ||, N pos = 6, DP = {lh, rh, lk, rk, la, ra} (71)

Mean Position and Orientation RMSE, Joint Angle RMSE and CC
In this experiment, multiple variations of the algorithm were tested as shown in Table 3. Firstly, L5S-3IMU is the algorithm described in this paper (Section 3) with parameters listed in Table 2. The parameter for L5S-3IMU were selected by taking the best joint CC (i.e., mean of free walk and dynamic movements) from a grid search of parameters σ 2 ω = {1, 10, 10 2 , 10 3 } rad 2 /s 2 and σ 2 ori = {10 −2 , 10 −1 , 1, 10} rad 2 . Secondly, CKF-3IMU and CKF-3IMU+D were the algorithms described in [15,31], respectively. Thirdly, CKF-3I-KB is a modified CKF-3IMU using similar parameters, measurement, and constraint functions as L5S-3IMU. The key difference between CKF-3IMU and CKF-3I-KB is that CKF-3I-KB allows knee bending, denoted by the suffix KB, during the constraint update. Fourthly, L5S-3I-NO is a variation of L5S-3IMU with σ 2 ω = 10 7 rad 2 /s 2 , σ 2 ori = 10 −1 rad 2 , and bω k = 0 rad. The parameters were chosen to have high uncertainty on the tracked orientation (i.e., effectively not using the orientation measurements at all), leading to a variation of L5S-3IMU that is similar to our prior work CKF-3IMU which assumed orientation measurements were noise-free. Lastly, the black box output (i.e., pelvis, thigh, and shank orientations) from the MVN Studio software (denoted as OSPS), which illustrates the performance of a widely-accepted commercial wearable IMC system with an OSPS configuration. For the first to fourth variations, the +D suffix means simulated inter-IMU distance measurements (σ dist = 0.1 m) was used instead of the pelvis height assumption. Table 3. The experiment was tested on the following algorithm variations.

Algorithm
Inter-IMU Distance Summary Description L5S-3IMU N Tracks position and orientation as described in Section 3 with parameters listed in Table 2. L5S-3IMU+D Y CKF-3IMU [15] N Only tracks position using a constrained KF. CKF-3IMU+D [31] Y CKF-3I-KB N Modified CKF-3IMU using similar parameters as L5S-3IMU (Table 2). This also allows knee bending during the constraint update. CKF-3I-KB+D Y L5S-3I-NO N L5S-3IMU with parameters that assume noise-free orientation (NO) measurements like CKF-3IMU. L5S-3I-NO+D Y OSPS N Output from a commercial OSPS wearable IMC system. Figure 4 shows the mean position and orientation RMSE, mean knee Y and hip joint angle RMSE (bias removed) and CC of different variations of CKF-3IMU and L5S-3IMU for both free walking and dynamic motions. Y, X, and Z refers to the sagittal, frontal, and transverse planes, respectively. CKF-3IMU performed well with free walking (e pos = 4.27 cm, e ori = 15.85 • , CC = 0.66) [15]. However, a more extensive evaluation showed that it performed poorly for certain dynamic movements (e.g., high-knee jog with e pos = 18.15 cm, e ori = 24.87 • , CC = 0.02). Removing the no-knee-bending assumption during the constraint update fixed this issue, as shown by the performance of CKF-3I-KB (e.g., high-knee jog improved by ∼9 cm e pos , ∼9 • e ori , ∼0.4 CC). L5S-3I-NO which is the L5S version of CKF-3IMU expectedly have similar performance with CKF-3I-KB (i.e., ∆e pos < 0.5 cm, ∆e ori < 1 • , and ∆CC 0.02 differences). L5S-3IMU, which tracked both position and orientation while assuming there is noise in the orientation measurements, had a slightly better performance (e.g., improved jumping jacks and high-knee jog by ∼0.1 CC, <0.03 CC difference with other movement types). The use of simulated distance measurement with σ dist = 0.1 m on CKF-3I-KB, L5S-3I-NO, and L5S-3IMU had slight effects for free walking, and a significant improvement for dynamic movements. For free walking, joint angle RMSE and CC of L5S-3IMU+D compared to L5S-3IMU improved by ∼1 • and <0.01 CC, while e pos and e ori slightly disimproved (<0.5 cm and <1 • ). The similar results suggest that inferring pelvis position from simulated distance measurement (σ dist = 0.1 m) is comparable to our pelvis height assumption at least for free walking. For dynamic movements, the e pos , e ori , joint angle RMSE, and CC of L5S-3IMU+D improved by 2-16 cm, 0-40 • , 1-9 • , and <0.42, respectively-more significantly for movements TUG and high-knee jog.  To give insight on how the accuracy of the simulated inter-IMU distance measurements affect pose estimation performance, Figure 5 shows the mean of knee Y and hip joint angle RMSE and CC at different σ dist values. At σ dist = 0.1 m, the simulation showed comparable performance between L5S-3IMU, which implements pelvis height assumption, and L5S-3IMU+D, which implements inter-IMU distance measurement to supplement the pelvis position estimate, for free walking. Significant improvement for dynamic movements can be seen even for σ dist = 0.2 m. These results suggest that the actual distance measurement sensor must have noise standard deviation σ dist ≤ 0.1 m to improve pose estimate performance. Note that the +D variation in Figure 4 and in the experiments that follow were evaluated at σ dist = 0.1 m. Figure 5. Joint angle RMSE (top) and CC (bottom) of free walk and dynamic movements at different noise level σ dist . The broken lines represent L5S-3IMU results (denoted as nD) where inter-IMU distance measurements were not used. The solid lines represent L5S-3IMU+D results (denoted as +D) where we can observe slight and great improvements for free walk and dynamic movements, respectively. Figure 6 shows the knee and hip joint angle RMSE (bias removed) and CC of L5S-3IMU and L5S-3IMU+D compared against the OMC output. Y, X, and Z refers to the sagittal, frontal, and transverse planes, respectively. Turning movements and half steps were manually removed from the per-step result of Walk movement and was denoted as Straight Walk. Note that sensor-to-body calibration was only done at the beginning of trial, not for each step. Between L5S-3IMU and L5S-3IMU+D, there was minimal hip and knee joint angle RMSE and CC improvement for free walking (∼1 • RMSE and ∼0.03 CC difference). However, there was significant improvement for most dynamic movements, specifically, speed-skater, jog, high-knee jog, and TUG (e.g., 4 • -17 • knee Y and hip Y joint angle RMSE improvements). Furthermore, the CC for dynamic movements started to reach similar performance with the free walk movement, indicating that inter-IMU distance measurements have indeed made the pose estimator capable of tracking more ADLs and not just walking.  Figure 7 shows a sample walk trial. At the peaks of knee Y angle, the distance between the pelvis and ankle positions of L5S-3IMU+D were a few cm shorter (i.e., pelvis position was lower than actual while ankle position was higher) than the actual distance resulting in higher knee Y angle peaks. Violations of our biomechanical constraints are also apparent at t = 4 to 5.5 s, where the subject makes a 180 • turn. After the turn, L5S-3IMU and L5S-3IMU+D were able to recover during the straight walking (t = 5.5 to 9.74 s of Figure 7). Notice that the bias between OSPS and OMC can be observed at t = 0 of the hip Y joint angle.  Table 4 shows the TTD, stride length, and gait speed accuracy computed from the global ankle position estimate of L5S-3IMU, L5S-3IMU+D, and the OMC system for free walk, jogging, and TUG. The use of inter-IMU distance measurements (σ dist = 0.1 m) helped improve the TTD, stride length, and gait speed accuracy of free walk and TUG (e.g., TTD improved from ∼9% to ∼5%). Refer to the code repository for links to videos of sample trials. where µ and σ denote mean and standard deviation. Error denotes estimate minus actual value, while TTD % dev denotes abs(error)/actual TTD.

Discussion
In this paper, a Lie group EKF algorithm for lower body pose estimation using only three IMUs, ergonomically placed on the ankles and sacrum to facilitate continuous recording outside the laboratory, was described and evaluated. The algorithm utilizes fewer sensors than other approaches reported in the literature, at the cost of reduced accuracy.

Mean Position and Orientation RMSE
The mean position and orientation RMSE of L5S-3IMU, L5S-3IMU+D, and related literature (sparse orientation poser (SOP) and sparse inertial poser (SIP) [17]) are listed in Table 5. SOP used orientation measured by IMUs and biomechanical constraints, while SIP used similar information but with the addition of acceleration. Both SOP and SIP were benchmarked against an OSPS system tracking the full body while our algorithm was benchmarked against an OMC system tracking only the lower body. The e pos and e ori (no bias) performance of L5S-3IMU and compared to SOP for free walking and jogging were comparable (∆e pos 0.1-0.5 cm and ∆e ori 2.5 • -3 • differences). The e pos and e ori (no bias) of SIP was better than L5S-3IMU and L5S-3IMU+D for free walking (∼2.1-2.5 cm and 6.5 • -7 • difference). Although this improvement was expected, as SIP optimizes the pose over multiple frames whereas our algorithm, like CKF-3IMU, optimizes the pose for each individual frame. For jumping jacks, the e ori of L5S-3IMU and L5S-3IMU+D was significantly (∼4 • -8 • ) better than SOP's and SIP's. However, this difference is probably because both SOP and SIP were evaluated on the full body (our algorithm was only evaluated on the lower body) and errors in arm pose estimation may have increased e ori for the SOP and SIP algorithms. Table 5. Mean position and orientation RMSE of L5S-3IMU, L5S-3IMU+D, OSPS, sparse orientation power (SOP) and sparse inertial poser (SIP) [17].  Comparing processing times, L5S-3IMU and L5S-3IMU+D were slower than CKF-3IMU, but can still be used in real-time; specifically, CKF-3IMU, L5S-3IMU, and L5S-3IMU+D processed a 1000-frame sequence (i.e., 10 s long) in ∼0.7, ∼2, ∼3.5 s, respectively, on an Intel Core i5-6500 3.2 GHz CPU [15], while SIP [17] took 7.5 min on a quad-core Intel Core i7 3.5 GHz CPU. All set-ups used single-core non-optimized Matlab code. Albeit slower than CKF-3IMU, L5S-3IMU and L5S-3IMU+D could also be used to provide real-time gait parameter measurement to inform actuation of assistive or rehabilitation robotic devices.

Hip and Knee Joint Angle RMSE and CC
The knee and hip joint angle RMSEs (no bias) and CCs of L5S-3IMU, L5S-3IMU+D, OSPS and related literature for straight walking (i.e., per step evaluation) are shown in Table 6 [7,15,37,38]. Similar to IMC based systems, L5S-3IMU and L5S-3IMU+D also follows the trend of having sagittal (Y axis) joint angles similar to that captured by OMC systems (0.95 knee Y and >0.83 hip Y CCs), but with significant difference in frontal and transverse (X and Z axis) joint angles [15,37]. CKF-3IMU performed slightly better (e.g., 0.03 knee Y, 0.09 hip Y CC), which is expected as the biomechanical constraint (i.e., no-knee-bending) assumption of CKF-3IMU was designed specifically for walking, at the cost of being less accurate for other dynamic movements. Both L5S-3IMU and L5S-3IMU+D were comparable, and at times even better (within 2.5 • RMSE, 0.1 CC difference) than the results of Hu et al. and Tadano et al., indicating excellent per-step reconstruction in the sagittal plane [7,38]. Hu et al. used 4 IMUs (two at the pelvis and one on each foot) and Tadano et al. used an OSPS configuration. Both systems can only estimate the pose in the sagittal plane.
Lastly, despite L5S-3IMU and L5S-3IMU+D achieving 0.95 joint angle CCs in the sagittal plane, the unbiased joint angle RMSE (>5 • ) makes its utility in clinical applications uncertain [41]. Although the algorithm is expected to work on pathological gait where our biomechanical assumptions are satisfied, overall performance still needs more improvement. To achieve clinical utility, one may either use more accurate sensors or average out cycle-to-cycle variation in estimation errors over many gait cycles; for example, use a more accurate distance measurement sensor (σ dist < 0.1 m). Furthermore, the accuracy must also be validated on a larger and more diverse cohort to quantify its ultimate clinical utility. The evaluation of how these solutions can bridge the gap to clinical application for the proposed system will be part of future work.

Spatiotemporal Gait Parameters
The focus of the proposed algorithms, L5S-3IMU and L5S-3IMU+D, are to estimate joint kinematics. However, as L5S-3IMU and L5S-3IMU+D both track the global position of the ankles, it is also capable of calculating spatiotemporal gait parameters (performance listed in Table 4). The TTD deviation of our algorithms compared against the gold standard OMC were not as good as CKF-3IMU [15] (3.6-3.81% TTD deviation) or other state-of-the-art dead reckoning algorithms [42,43] (0.2-1.5% TTD deviation). Two possible sources of inaccuracy lies (1) in the dead reckoning approximation done in the prediction step, and (2) in the assumption that the velocity of the shank IMU is zero when the associated foot touches the floor, but of course this IMU continues to move with some small velocity on the lower shank during the stance phase. To illustrate the dead reckoning approximation, let us look at the predicted pelvis pose in Equation (73). In our algorithm, we assumed ψ p ≈ I 3×3 (note that Φ(−∆t pω k ) ≈ I 3×3 andR p+ k−1 (Ȓ p k ) T ≈ I 3×3 since ∆t pω k is small) which did not significantly affect the joint kinematic estimate, but slightly affected the global position estimate. Nevertheless, body drift has been reduced substantially compared to Marcard et al.'s SIP [17].

Limitations and Future Work
L5S has similar pelvis drift, covariance matrix numerical issue, and flat floor limitation as CKF-3IMU, which is expected as L5S implements the same measurement and constraint update as CKF-3IMU, albeit formulated using Lie group representation instead of vectors and quaternions [15]. The pelvis height and flat floor assumption helps prevent the pelvis and the ankles from drifting towards each other (i.e., pelvis drift downward while ankles drift upward). However, it will also prevent accurate pose estimation of motions such as sitting, lying down, or standing on one leg, where the pose is maintained for a duration much longer than that of a typical gait cycle. The covariance limiter (Section 3.2.2) helps prevent the covariance becoming badly conditioned (i.e., singular), especially for longer duration trials (e.g., 5-minute walk) where the position uncertainty grows at a faster rate for the pelvis position than the ankle position. As can be observed from Figure 6, substituting the pelvis height assumption with inter-IMU distance measurements can increase the algorithm's accuracy especially for tracking dynamic movements. If the distance measurement is accurate enough (i.e., smaller σ 2 dist ), the inter-IMU distance measurement update may be enough to limit the growth of pelvis position uncertainly and possibly making the covariance limiter not needed. Figure 6 shows that the optimized performance of L5S-3IMU, even if it allows the tracked orientation to be corrected by inter-IMU distance measurements and the tracked position estimate, was only slightly better than CKF-3IMU/L5S-3I-NO, which effectively assumed the measurement input from the orientation estimation algorithm to be perfect (i.e., trusted the tracked orientation less). As L5S-3IMU requires more computing resources, such result suggests that CKF-3IMU may be more suitable to use when computing power is limited. To fully leverage the advantages brought by the Lie group representation, additional sensor measurements that can help correct tracked orientation will be needed (e.g., estimating angle of arrival between two sensors [44] or using fish eye cameras to improve pose estimate [45]).
Additional sensor measurements provide new opportunities for automatic calibration even under RSC configuration. IMC systems typically need anthropometric measurements (i.e., measurement of body segments such as d ls ) beforehand. By taking the initial distance measurement at some predetermined posture, anthropometric measurements can be automatically inferred. The formulation for a hinge joint with two IMUs on both sides has been leveraged to enable automatic sensor-to-segment calibration (i.e., align sensor frame to body frame) and even a completely magnetometer free orientation estimation [46,47]. Magnetometer free orientation estimation rids us of the yaw offset issue from an inhomogeneous magnetic field in indoor environments, typically with stronger disturbances closer to the floor [48]. An approach using a hinge joint with two IMUs may not be applicable to RSC configurations (e.g., our algorithm only has one IMU on one side of the hinge joint). However, distance measurements may be use to compensate for the missing IMU information from the uninstrumented segment, and a modified version may be developed for a RSC configuration.
Enabling longer-term tracking of ADL in the subject's natural environment may lead to novel investigations of movement disorder progression and the identification of early intervention opportunities. This work is just one of the early steps towards seamless remote gait monitoring.
Developing solutions to further increase accuracy, increase the number of body segments tracked (e.g., track full body under RSC [17]), or use even fewer IMUs (tracking lower body using two IMUs [49]) will be investigated in the future.

Conclusions
This paper presented a Lie group CEKF-based algorithm (L5S-3IMU) to estimate lower limb kinematics using a RSC configuration of IMUs, supplemented by inter-IMU distance measurements in one implementation. The knee and hip joint angle RMSEs in the sagittal plane for straight walking were 7.6 • ± 2.6 • and 6.6 • ± 2.7 • , respectively, while the CCs were 0.95 ± 0.03 and 0.87 ± 0.16, respectively. We also showed that inter-IMU distance measurement is a promising new source of information to improve the pose estimation of IMC under a RSC configuration. Simulations show that performance improved dramatically for dynamic movements even at higher noise levels (e.g., σ dist = 0.2 m), and that similar performance to L5S-3IMU was achieved at σ dist = 0.1 m for free walk movements. However, further validation is recommended with actual distance measurement from real sensors. The source code for the L5S algorithm, and links to sample videos will be made available at https://git.io/JTRQ3. First, we solve for an estimated left knee angle,θ lk k (Equation (47)), from the measured pelvis to left ankle distance,d pla k . The pelvis to left ankle vector, τ pla m (μ − k , θ lk k ) (Equation (A6)), can be defined as the sum of the mid-pelvis to hip, thigh long axis, and shank long axis vectors. τ pla (μ − k , θ lk k ) = ψ pla =half pelvis y-axis + shank z-axis d p 2T pk i y − d lsT lsk i z +d ltT lsk thigh z-axis in shank frame (i x sin (θ lk k ) − i z cos (θ lk k )) (A1) By definition of (d Solving forθ lk k from Equation (A3) gives us a quadratic equation with two solutions as shown in Equations (A5) and (47). Between the two solutions,θ lk k is set as theθ lk k whose value is closer to the current left knee angle estimate from the prediction step. This solution serves as a pseudomeasurement of the knee angle.θ lk k = cos −1 αγ±β √ α 2 +β 2 −γ 2 α 2 +β 2 (A5) Finally, Z pla,k , the KF measurement shown in Eqs. (A6) and (46), is the inter-IMU vector between the pelvis and left ankle, calculated using Equation (A1) with inputθ lk k .