Contact-Aided Invariant Extended Kalman Filtering for Robot State Estimation

TL;DR

This paper introduces a contact-aided invariant extended Kalman filter for legged robot state estimation, leveraging Lie group theory to improve convergence, accuracy, and robustness over traditional methods.

Contribution

It develops a novel InEKF that exploits system symmetries, enabling better convergence and accuracy in robot pose and velocity estimation.

Findings

Outperforms quaternion-based EKF in simulations and experiments

Demonstrates improved convergence and robustness

Effectively incorporates contact and inertial data

Abstract

Legged robots require knowledge of pose and velocity in order to maintain stability and execute walking paths. Current solutions either rely on vision data, which is susceptible to environmental and lighting conditions, or fusion of kinematic and contact data with measurements from an inertial measurement unit (IMU). In this work, we develop a contact-aided invariant extended Kalman filter (InEKF) using the theory of Lie groups and invariant observer design. This filter combines contact-inertial dynamics with forward kinematic corrections to estimate pose and velocity along with all current contact points. We show that the error dynamics follows a log-linear autonomous differential equation with several important consequences: (a) the observable state variables can be rendered convergent with a domain of attraction that is independent of the system's trajectory; (b) unlike the standard EKF, neither the linearized error dynamics nor the linearized observation model depend on the current state estimate, which (c) leads to improved convergence properties and (d) a local observability matrix that is consistent with the underlying nonlinear system. Furthermore, we demonstrate how to include IMU biases, add/remove contacts, and formulate both world-centric and robo-centric versions. We compare the convergence of the proposed InEKF with the commonly used quaternion-based EKF though both simulations and experiments on a Cassie-series bipedal robot. Filter accuracy is analyzed using motion capture, while a LiDAR mapping experiment provides a practical use case. Overall, the developed contact-aided InEKF provides better performance in comparison with the quaternion-based EKF as a result of exploiting symmetries present in system.

Tables

Table 1. Table 1 : Experiment Discrete Noise Statistics and Initial Covariance

Measurement Type noise st. dev. Linear Acceleration $0.04m/{\sec }^2$ Angular Velocity $0.002rad/\sec$ Accelerometer Bias $0.001m/{\sec }^3$ Gyroscope Bias $0.001rad/{\sec }^2$ Contact Linear Velocity $0.05m/\sec$ Joint Encoders $1.0\mathrm{deg}$

State Element initial st. dev. Orientation of IMU $30.0\mathrm{deg}$ Velocity of IMU $1.0m/\sec$ Position of IMU $0.1m$ Position of Right Foot $0.1m$ Position of Left Foot $0.1m$ Gyroscope Bias $0.005rad/\sec$ Accelerometer Bias $0.05m/{\sec }^2$

Table 2. Table 2 : Summary of World-centric State Estimator

State Definition	Deterministic Nonlinear Dynamics

Equations

$${(\cdot )}^{\wedge }:{\mathbb{R}}^{\dim \mathfrak{g}}\to \mathfrak{g}$$

$$\exp (\mathbf{\xi })={\exp }_m({\mathbf{\xi }}^{\wedge }),$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\text{X}}_t=f_{u_t}({\text{X}}_t),$$

$$\begin{array}{rl}{\mathbf{\eta }}_t^r& ={\bar{\text{X}}}_t{\text{X}}_t^{-1}=({\bar{\text{X}}}_t\text{L})({\text{X}}_t\text{L}{)}^{-1}\text{(Right-Invariant)}\\ {\mathbf{\eta }}_t^l& ={\text{X}}_t^{-1}{\bar{\text{X}}}_t=(\text{L}{\bar{\text{X}}}_t{)}^{-1}(\text{L}{\text{X}}_t),\text{(Left-Invariant)}\end{array}$$

$$f_{u_t}({\text{X}}_1{\text{X}}_2)=f_{u_t}({\text{X}}_1){\text{X}}_2+{\text{X}}_1{f}_{u_t}({\text{X}}_2)-{\text{X}}_1{f}_{u_t}({\text{I}}_d){\text{X}}_2$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\eta }}_t^r$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\eta }}_t^l$$

$$g_{u_t}(\exp (\mathbf{\xi }))\triangleq {({\text{A}}_t\mathbf{\xi })}^{\wedge }+\mathcal{O}(||\mathbf{\xi }|{|}^2).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\xi }}_t={\text{A}}_t{\mathbf{\xi }}_t.$$

$${\mathbf{\eta }}_t=\exp ({\mathbf{\xi }}_t);$$

$${\text{Y}}_t$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\left[\begin{array}{c}{q}_x\\ q_y\\ q_z\end{array}\right]=\left[\begin{array}{ccc}1& \sin (q_x)\tan (q_y)& \cos (q_x)\tan (q_y)\\ 0& \cos (q_x)& -\sin (q_x)\\ 0& \sin (q_x)\sec (q_y)& \cos (q_x)\sec (q_y)\end{array}\right]\left[\begin{array}{c}{w}_x\\ w_y\\ w_z\end{array}\right].$$

$$\frac{\mathrm{d}}{\mathrm{d}t}\delta {\text{q}}_t=g(\delta {\text{q}}_t,{\tilde{\mathbf{\omega }}}_t,{\text{q}}_t).$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}\delta {\text{q}}_t& \approx \left[\begin{array}{ccc}0& ({\omega }_z{\bar{c}}_x+{\omega }_y{\bar{s}}_x)/{\bar{c}}_y^2& {\bar{t}}_y({\omega }_y{\bar{c}}_x-{\omega }_z{\bar{s}}_x)\\ 0& 0& {\omega }_z{\bar{c}}_x+{\omega }_y{\bar{s}}_x\\ 0& ({\bar{s}}_y({\omega }_z{\bar{c}}_x+{\omega }_y{\bar{s}}_x))/{\bar{c}}_y^2& ({\omega }_y{\bar{c}}_x-{\omega }_z{\bar{s}}_x)/{\bar{c}}_y\end{array}\right]\delta {\text{q}}_t\\ & \triangleq \text{A}({\tilde{\mathbf{\omega }}}_t,{\bar{\text{q}}}_t)\delta {\text{q}}_t,\end{array}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\text{R}}_t={\text{R}}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times },$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\eta }}_t& ={\text{R}}_t^{\mathsf{T}}{\bar{\text{R}}}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }+{\left({\text{R}}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }\right)}^{\mathsf{T}}{\bar{\text{R}}}_t\\ & ={\text{R}}_t^{\mathsf{T}}{\bar{\text{R}}}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\text{R}}_t^{\mathsf{T}}{\bar{\text{R}}}_t\\ & ={\mathbf{\eta }}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\mathbf{\eta }}_t\\ & =g({\mathbf{\eta }}_t,{\tilde{\mathbf{\omega }}}_t).\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}(\exp ({\mathbf{\xi }}_t))& =\exp ({\mathbf{\xi }}_t){\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }\exp ({\mathbf{\xi }}_t)\\ \frac{\mathrm{d}}{\mathrm{d}t}\left(\text{I}+{\left({\mathbf{\xi }}_t\right)}_{\times }\right)& \approx \left(\text{I}+{\left({\mathbf{\xi }}_t\right)}_{\times }\right){\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }\left(\text{I}+{\left({\mathbf{\xi }}_t\right)}_{\times }\right)\\ ⟹\frac{\mathrm{d}}{\mathrm{d}t}{\left({\mathbf{\xi }}_t\right)}_{\times }& ={\left({\mathbf{\xi }}_t\right)}_{\times }{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\left({\mathbf{\xi }}_t\right)}_{\times }\\ & ={\left(-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\mathbf{\xi }}_t\right)}_{\times }\\ ⟹\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\xi }}_t& ={\left(-{\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\mathbf{\xi }}_t\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\eta }}_t& ={\text{R}}_t^{\mathsf{T}}{\mathbf{\eta }}_0{\text{R}}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\text{R}}_t^{\mathsf{T}}{\mathbf{\eta }}_0{\text{R}}_t\\ & ={\mathbf{\eta }}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\mathbf{\eta }}_t\end{array}$$

$$\begin{array}{rl}{\mathbf{\eta }}_t& ={\text{R}}_t^{\mathsf{T}}{\mathbf{\eta }}_0{\text{R}}_t\\ ⟹\exp ({\mathbf{\xi }}_t)& ={\text{R}}_t^{\mathsf{T}}\exp ({\mathbf{\xi }}_0){\text{R}}_t=\exp ({\text{R}}_t^{\mathsf{T}}{\mathbf{\xi }}_0)\\ ⟹{\mathbf{\xi }}_t& ={\text{R}}_t^{\mathsf{T}}{\mathbf{\xi }}_0\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\xi }}_t& =-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\text{R}}_t^{\mathsf{T}}{\mathbf{\xi }}_0\\ & =-{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }{\mathbf{\xi }}_t\end{array}$$

$${\text{X}}_t\triangleq \left[\begin{array}{cccccc}{\text{R}}_{\text{WB}}(t)& {}_{\text{W}}{\text{v}}_{\text{B}}(t)& {}_{\text{W}}{\text{p}}_{\text{WB}}(t)& {}_{\text{W}}{\text{p}}_{{\text{WC}}_1}(t)& \cdots & {}_{\text{W}}{\text{p}}_{{\text{WC}}_N}(t)\\ {\text{0}}_{1,3}& 1& 0& 0& \cdots & 0\\ {\text{0}}_{1,3}& 0& 1& 0& \cdots & 0\\ {\text{0}}_{1,3}& 0& 0& 1& \cdots & 0\\ ⋮& ⋮& ⋮& ⋮& \ddots & ⋮\\ {\text{0}}_{1,3}& 0& 0& 0& \cdots & 1\end{array}\right].$$

$$\begin{array}{r}{\text{X}}_t\triangleq \left[\begin{array}{cccc}{\text{R}}_t& {\text{v}}_t& {\text{p}}_t& {\text{d}}_t\\ {\text{0}}_{1,3}& 1& 0& 0\\ {\text{0}}_{1,3}& 0& 1& 0\\ {\text{0}}_{1,3}& 0& 0& 1\end{array}\right],{\text{u}}_t=\left[\begin{array}{c}{}_{\text{B}}{\tilde{\mathbf{\omega }}}_{\text{WB}}(t)\\ {}_{\text{B}}{\tilde{\text{a}}}_{\text{WB}}(t)\end{array}\right]\triangleq \left[\begin{array}{c}{\tilde{\mathbf{\omega }}}_t\\ {\tilde{\text{a}}}_t\end{array}\right],\end{array}$$

$${\mathbf{\xi }}^{\wedge }=\left[\begin{array}{cccc}{\left({\mathbf{\xi }}^R\right)}_{\times }& {\mathbf{\xi }}^v& {\mathbf{\xi }}^p& {\mathbf{\xi }}^d\\ {\text{0}}_{1,3}& 0& 0& 0\\ {\text{0}}_{1,3}& 0& 0& 0\\ {\text{0}}_{1,3}& 0& 0& 0\end{array}\right],$$

$${\mathrm{Ad}}_{{\text{X}}_t}=\left[\begin{array}{cccc}\text{R}& \text{0}& \text{0}& \text{0}\\ {\left({\text{v}}_t\right)}_{\times }{\text{R}}_t& {\text{R}}_t& \text{0}& \text{0}\\ {\left({\text{p}}_t\right)}_{\times }{\text{R}}_t& \text{0}& {\text{R}}_t& \text{0}\\ {\left({\text{d}}_t\right)}_{\times }{\text{R}}_t& \text{0}& \text{0}& {\text{R}}_t\end{array}\right].$$

$${\tilde{\mathbf{\omega }}}_t$$

$${\tilde{\text{a}}}_t$$

$${}_{\text{W}}{\tilde{\text{v}}}_{\text{C}}={\text{0}}_{3,1}={}_{\text{C}}{\text{v}}_{\text{C}}+{\text{w}}_t^v,{\text{w}}_t^v\sim \mathcal{G}\mathcal{P}\left({\text{0}}_{3,1},{\mathbf{\Sigma }}^v\delta (t-t^{\prime })\right).$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\text{R}}_t& ={\text{R}}_t{\left({\tilde{\mathbf{\omega }}}_t-{\text{w}}_t^g\right)}_{\times }\\ \frac{\mathrm{d}}{\mathrm{d}t}{\text{v}}_t& ={\text{R}}_t({\tilde{\text{a}}}_t-{\text{w}}_t^a)+\text{g}\\ \frac{\mathrm{d}}{\mathrm{d}t}{\text{p}}_t& ={\text{v}}_t\\ \frac{\mathrm{d}}{\mathrm{d}t}{\text{d}}_t& ={\text{R}}_t{\text{h}}_R({\tilde{\mathbf{\alpha }}}_t)(-{\text{w}}_t^v),\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\text{X}}_t& =\left[\begin{array}{cccc}{\text{R}}_t{\left({\tilde{\mathbf{\omega }}}_t\right)}_{\times }& {\text{R}}_t{\tilde{\text{a}}}_t+\text{g}& {\text{v}}_t& {\text{0}}_{3,1}\\ {\text{0}}_{1,3}& 0& 0& 0\\ {\text{0}}_{1,3}& 0& 0& 0\\ {\text{0}}_{1,3}& 0& 0& 0\end{array}\right]-\left[\begin{array}{cccc}{\text{R}}_t& {\text{v}}_t& {\text{p}}_t& {\text{d}}_t\\ {\text{0}}_{1,3}& 1& 0& 0\\ {\text{0}}_{1,3}& 0& 1& 0\\ {\text{0}}_{1,3}& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}{\left({\text{w}}_t^g\right)}_{\times }& {\text{w}}_t^a& {\text{0}}_{3,1}& {\text{h}}_R({\tilde{\mathbf{\alpha }}}_t){\text{w}}_t^v\\ {\text{0}}_{1,3}& 0& 0& 0\\ {\text{0}}_{1,3}& 0& 0& 0\\ {\text{0}}_{1,3}& 0& 0& 0\end{array}\right]\\ & \triangleq f_{u_t}({\text{X}}_t)-{\text{X}}_t{{\text{w}}_t}^{\wedge },\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\eta }}_t^r& =f_{u_t}({\mathbf{\eta }}_t^r)-{\mathbf{\eta }}_t^r{f}_{u_t}({\text{I}}_d)+({\bar{\text{X}}}_t{{\text{w}}_t}^{\wedge }{\bar{\text{X}}}_t^{-1}){\mathbf{\eta }}_t^r\\ & \triangleq g_{u_t}({\mathbf{\eta }}_t^r)+{{\bar{\text{w}}}_t}^{\wedge }{\mathbf{\eta }}_t^r\end{array}$$

$$\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\xi }}_t& ={\text{A}}_t^r{\mathbf{\xi }}_t+{\bar{\text{w}}}_t={\text{A}}_t{\mathbf{\xi }}_t+{\mathrm{Ad}}_{{\bar{\text{X}}}_t}{\text{w}}_t\\ {\mathbf{\eta }}_t^r& =\exp ({\mathbf{\xi }}_t).\end{array}$$

Full text

1 Introduction

Legged robots have the potential to transform the logistics and package delivery industries, become assistants in our homes, and aide in search and rescue [29]. In order to develop motion planning algorithms and robust feedback controllers for these tasks, accurate estimates of the robot’s state are needed. Some states, such as joint angles, can be directly measured using encoders, while other states, such as the robot’s pose and velocity, require additional sensors. Most legged robots are equipped with an inertial measurement unit (IMU) that can measure linear acceleration and angular velocity, albeit with noise and bias perturbations. Consequently, nonlinear observers are typically used to fuse leg odometry and inertial measurements to infer trajectory, velocity, and calibration parameters [68, 14, 35, 53]. In view of a practical solution, designing a globally convergent observer is sacrificed for one with at best local properties, such as the extended Kalman filter (EKF) [42, 52, 79, 75]. This EKF-based approach is computationally efficient and easily customizable, allowing successful implementation on a number of legged robots with rigorous real-time performance requirements [15, 17, 36, 65].

Accurate pose estimation can be combined with visual data to build maps of the environment [37]. Then such maps can be used in gait selection to improve the stability of a robot while walking on uneven terrains and as a basis for high-level motion planning. Although there have been many recent advancements in visual-inertial-odometry and simultaneous localization and mapping (SLAM) [66, 48, 38], these algorithms often rely on visual data for pose estimation. This means that the observer (and ultimately the feedback controller) can be adversely affected by rapid changes in lighting as well as the operating environment. It is therefore beneficial to develop a low-level state estimator that fuses data only from proprioceptive sensors to form accurate high-frequency state estimates. This approach was taken by Bloesch et al. [15] when developing a quaternion-based extended Kalman filter (QEKF) that combines inertial, contact, and kinematic data to estimate the robot’s base pose, velocity, and a number of contact states. In this article, we expand upon these ideas to develop an invariant extended Kalman filter (InEKF) that has improved convergence and consistency properties allowing for a more robust observer that is suitable for long-term autonomy.

The theory of invariant observer design is based on the estimation error being invariant under the action of a matrix Lie group [1, 20], which has recently led to the development of the InEKF111We use the InEKF acronym to distinguish from an iterated-EKF (IEKF). [18, 8, 10, 11] with successful applications and promising results in simultaneous localization and mapping [8, 86] and aided inertial navigation systems [4, 5, 8, 83]. The invariance of the estimation error with respect to a Lie group action is referred to as a symmetry of the system [5]. Summarized briefly, Barrau and Bonnabel [10] showed that if the state is defined on a Lie group, and the dynamics satisfy a particular “group affine” property, then the symmetry leads to the estimation error satisfying a “log-linear” autonomous differential equation on the Lie algebra. In the deterministic case, this linear system can be used to exactly recover the estimated state of the nonlinear system as it evolves on the group. The log-linear property therefore allows the design of a nonlinear observer or state estimator with strong convergence properties.

1.1 Contribution

In this article, we derive an InEKF for a system containing IMU and contact sensor dynamics, with forward kinematic correction measurements. We show that the defined deterministic system satisfies the “group affine” property, and therefore, the error dynamics is exactly log-linear. In practice, with the addition of sensor noise and IMU bias, this log-linear error system is only approximate; we show, however, that in many cases the proposed InEKF is still preferred over standard QEKFs due to superior convergence and consistency properties. We demonstrate the utility and accuracy of the developed observer through a series of LiDAR mapping experiments performed on a Cassie-series biped robot.

In summary, this work makes the following contributions:

Derivation of a continuous-time right-invariant EKF for an IMU/contact process model with a forward kinematic measurement model; the observability analysis is also presented; 2. 2.

State augmentation with IMU biases; 3. 3.

Evaluations of the derived observer in simulation and hardware experiments using a 3D bipedal robot; 4. 4.

Alternative derivation of the observer using a left-invariant error definition; 5. 5.

Detailed explanation of the connection between the invariant error choice and the world-centric and robo-centric estimator formulations; 6. 6.

Equations provided for analytical discretization of the proposed observers; and 7. 7.

Development of an open-source C++ library for aided-inertial navigation using the InEKF https://github.com/RossHartley/invariant-ekf.

1.2 Outline

The remainder of this article is organized as follows. Background and related work are given in section 2. Section 3 provides the necessary preliminary material needed for understanding the InEKF formulation, which is motived by an example from attitude dynamics in section 4. Section 5 provides the derivation of a right-invariant extended Kalman filter (RIEKF) for contact-inertial navigation with a right-invariant forward kinematic measurement model. In Section 6, we present simulation results comparing the convergence properties to a state-of-the-art QEKF. Section 7 discusses the state augmentation of the previously derived InEKF with IMU bias. The consequences of adding and removing of contact points in the estimator are described in Section 8. Experimental evaluations on a 3D biped robot, shown in Figure 1, are presented in Section 9 along with an application towards LiDAR-based terrain mapping. Section 10 provides an alternative derivation of the proposed observer using the left-invariant dynamics along with an explanation of how to easily switch between the two formulations. Section 11 details how these equations can be modified to create a “robot-centric” estimator. Section 12 itemizes additional sensor measurements that fit within the InEKF framework and describes the relation between the developed filter and landmark-based SLAM. Finally, Section 13 concludes the article and suggests future directions. Details about time-discretization and useful Lie group expressions are given in the appendix.

2 Background and Related work

In this section, we first review the Kalman filtering literature to locate the proposed state estimator within the relate work. Then, we review the state estimation techniques for humanoid and legged robots.

2.1 Kalman Filtering

Filtering methods involve estimating the robot’s current state (and potentially landmarks) using the set of all measurements up to the current time [2, 41, 73]. When the process model and measurements are linear, and the noise is white and Gaussian, Kalman filtering [51] provides an optimal method (minimum mean squared error) for state estimation. The general process for Kalman filtering involves two phases, propagation and correction. The state is typically represented using a Gaussian random vector, which is parameterized by a mean and a covariance. During the propagation phase, the previous state and covariance estimate are propagated forward in time using the system dynamics (alternatively known as the process model). When a measurement is obtained, the state and covariance estimate are corrected using the measurement model along with an associated measurement noise covariance.

Although the Kalman filter provides a method for optimal linear filtering, most practical mobile robots have nonlinear system dynamics, and many useful sensor models are also nonlinear. For these cases, an EKF can be designed, which utilizes Taylor series expansions to linearize the process and measurement models around the current state estimate [77]. Due to its low computational complexity and accurate performance, the EKF quickly became the de facto standard of nonlinear filtering for many mobile robot applications, including wheeled vehicles, drones, and legged robots [24, 70, 15]. The EKFs has also been proposed to solve the SLAM problem [74]. However, since the nonlinear system is linearized about the current state estimate, the EKF is, at best, only a locally stable observer [76, 52]. The local convergence proofs are based on Lipschitz bounds of the nonlinear terms in a model, and hence “the more nonlinear a system is, the worse an EKF may perform”. Importantly, if the state estimate is initialized poorly, it is possible for the filter to diverge. In addition, because an EKFs uses a systems linearization about the current estimate, states that are unobservable can spuriously be treated as observable by the filter. While, this can be mitigated through the use of an observability-constrained EKF developed by Huang et al. [49], it cannot be altogether avoided.

For many systems, the error-state (or indirect) EKF offers superior performance to the standard (total state or direct) form. As the name implies, the error-state extended Kalman filter (ErEKF) is formulated using the errors, such as pose and velocity errors, as the filter variables, while the standard EKF tracks the states themselves (pose and velocity directly) [59, 69, 79, 75]. Under small noise assumptions, this leads to linear error dynamics which are then used for covariance propagation in the error-state filter. The measurement model is also rewritten with respect to these errors. Although the error dynamics are linear in the error variables, they may still depend on the current state estimate, which if initialized poorly, will degrade the performance of the filter. However, the approximately linear nature of the error dynamics may respect the linear assumptions of the original Kalman filter better than the underlying system dynamics, which can lead to improved performance [59].

Perhaps the most important feature of the ErEKF is the ability to circumvent dynamic modeling [69]. This is done by replacing the potentially complicated process model with a relatively simple IMU integration model (also known as strapdown modeling) [78, 62, 82]. Essentially, the IMU’s angular velocity and linear acceleration measurements are integrated to propagate the state estimate, while the covariance is propagated using the error dynamics. Additional (independent) sensor measurements will correct the estimated error, which can then be used to update the state estimate. Using this method, there is no longer a need to formulate complicated, platform-specific dynamics models, which may require a large number of state variables and is likely to be exceedingly nonlinear. The “strapdown” ErEKFs has proven to yield highly accurate results (even with a low-cost IMUs) and continues to form the basis of many inertial navigation systems [13, 54, 79, 75, 15].

In the standard formulation of Kalman filtering theory, the system evolves on Euclidean spaces. However, in many cases, the state variables we are interested in lie on a manifold. For example, the orientation of a 3D rigid body is represented by an element of the special orthogonal group, $\mathrm{SO}(3)$. This matrix Lie group is defined by the set of orthogonal $3\times 3$ matrices with a determinant of one. While the matrix contains nine variables, yet the dimension of the manifold is only three. One common approach is to parameterize $\mathrm{SO}(3)$ using local coordinates such as three Euler or Tait-Bryan angles [39, 71]. This allows the standard ErEKF equations to be applied; however, these local parameterizations are often arbitrary (and therefore confusing) and contain singularities (the well-known Gimbal lock problem). Alternatively, it is possible to represent 3D orientation using quaternions, which are a four-dimensional double cover (a two to one diffeomorphism) of $\mathrm{SO}(3)$. Using quaternions eliminates the singularities; however, modifications to the standard ErEKF equations have to be made [75, 79]. In brief, while the orientation is represented by a four-dimensional quaternion, the orientation error has to be defined by a $3-$vector (in the Lie algebra of $\mathrm{SO}(3)$) and the associated covariance by a $3\times 3$ matrix in order to prevent degeneracy. Also, the orientation corrections are done through quaternion multiplication instead of vector addition. This QEKF 222The quaternion-based formulation of EKFs is also sometimes called “multiplicative filtering” (MEKF) due to the orientation correction being done through quaternion multiplication [60]. has been well studied and implemented on a number of platforms, ranging from spacecraft [63, 54] to legged robots [15, 68, 36].

It turns out, many useful robot states can be characterized using matrix Lie groups. Examples include 2D orientation, $\mathrm{SO}(2)$, 3D orientation, $\mathrm{SO}(3)$, and 3D pose (orientation and position), $\mathrm{SE}(3)$. If the state to be estimated is a matrix Lie group, it is possible to further improve the EKF filtering approach. Bourmaud et al. [22, 23] developed versions of both discrete and continuous-time EKFs for systems where the state dynamics and measurements evolve on matrix Lie groups. In these formulations, noise is represented as a concentrated Gaussian on Lie groups [80, 81], which is a generalization of the multivariate Gaussian distribution. In essence, noise is represented as a Gaussian in the tangent space about a point on the manifold. This noise is then mapped to the Lie group through the use of the group’s exponential map, resulting in a decidedly non-Gaussian distribution on the manifold. An improved state estimate is obtained due to the filter taking into account the geometry and structure of the problem [23].

Most recently, a new type of EKF has been developed that is rooted in the theory of invariant observer design, in which the estimation error is invariant under the action of a Lie group [1, 20]. This invariance is referred to as the symmetries of the system [5]. This work led to the development of the InEKF [18, 8, 10, 11], with successful applications and promising results in SLAM [8, 86] and aided INSs [21, 4, 5, 8, 83]. Similar to the above mentioned EKFs on matrix Lie groups, the state is again represented as a matrix Lie group and the noise as a concentrated Gaussian on the group.

However, the InEKF exploits available system symmetries to further improve filtering results. The culminating result of the InEKF states that if a system satisfies a “group-affine” property, the estimation error satisfies a “log-linear” autonomous differential equation on the Lie algebra of the corresponding Lie group [11, 8]. In other words, the system linearization does not depend upon the estimated states. Therefore, one can design a nonlinear state estimator with strong convergence properties. Surprisingly, many mobile robot state estimation problems can be solved within the InEKF framework. This includes attitude estimation [18], inertial odometry [11, 8], velocity-aided inertial navigation [21], landmark-aided navigation [8], GPS and magnetometer-aided navigation [4], and even EKF-based SLAM. In this article, we extend this class of solutions to contact-aided inertial navigation [46], where forward kinematics is used to correct an inertial and contact-based process model. This approach successfully allows an InEKF to be used for legged robot state estimation.

2.2 Legged Robot and Humanoid State Estimation

Legged robots are a subclass of mobile robots that locomote through direct and switching contact with the environment. These robots typically contain proprioceptive sensors, such as IMU, joint encoders, and contact sensors. In addition, some legged robots, especially humanoids, also have access to exteroceptive sensors, namely cameras and LiDARs. As with all mobile robots, state estimation for legged robots is critical for mapping, planning, designing feedback controllers, and developing general autonomy. In this section, an overview of notable techniques for legged robot state estimation is given.

The simplest approach for estimation of a legged robot’s spatial location and velocity is kinematic dead-reckoning, otherwise known as kinematic odometry. This involves estimating relative transformations using only kinematic and contact measurements. In particular, encoder measurements and the kinematics model are used to track the position, orientation, and velocity of the robot’s base frame based on the assumption that a stance foot remains fixed to the ground. Although this method can be easily implemented, the state estimate is typically noisy due to kinematic modeling errors, encoder noise, and foot slip [67]. When only one foot is in contact with the ground, this “static contact assumption” may be violated. For example, if the robot has point feet, the stance foot position may remain fixed, but the foot orientation is free to rotate (without changing the joint angles). Therefore, a gyroscope is often used to provide angular velocity measurements which allows the robot’s body velocity to be recovered.

Alternatively, if the terrain is known a priori and at least three noncollinear point feet are on the ground, Lin et al. [55] showed that the robot’s instantaneous base pose can be computed through kinematics. These kinematic-based methods have been implemented on a number of legged robots including a planar one-legged hopper [47], the CMU Ambler hexapod [67], the RHex hexapod [55], and the biped robot MARLO [31, 32]. However, due to the high amounts of noise coming from encoders and foot slip, the velocity estimate typically needs to be heavily filtered before becoming usable in the feedback controller [44]. In addition, this noise causes the position and orientation estimates to drift substantially rendering the estimator impractical for mapping and autonomy tasks.

Fortunately, legged robots are often equipped with additional sensors such as IMUs, GPS, cameras, or LiDARs which provide independent, noisy odometry measurements. Much of the literature on legged robot state estimation focuses on fusing these measurements (potentially with kinematic odometry) using filtering and smoothing methods. Singh et al. [72] combined inertial measurements with optical flow measurements in a four-phase hybrid EKF. This required explicit dynamic modeling of the robot in flight, landing, stance, and thrust phases. Lin et al. [56] took a similar model-based approach and used an EKF to fuse kinematic information with IMU measurements to estimate the state of a hexapod. Cobano et al. [30] developed an EKF that fuses kinematic odometry and magnetometer readings with position measurements from a GPS to localize a SILO4 quadruped outdoors. This implementation fixes the issues with unbounded drift, but cannot operate in GPS-denied environments. If a prior terrain map is known, Chitta et al. [28] showed that it is possible to solve the localization problem for legged robots using only proprioceptive sensors and a particle filter. The key idea was that if the robot “senses” that a terrain change through kinematics, then this limits the potential locations the robot can be in a known map. The method was demonstrated on the LittleDog quadruped.

A breakthrough came in 2012 when Bloesch et al. [15] combined inertial and kinematic measurements in an observability-constrained ErEKF using the strapdown IMU modeling approach. In this work, no a priori knowledge of the terrain is assumed, and the IMU integration model completely eliminates the need for dynamic modeling of the robot. Therefore, the derived filter equations are general enough to be used on any legged robot. The key idea was to augment the state vector with the set of all foot positions currently in contact with the environment. During the prediction phase, the foot contact dynamics are assumed to be Brownian motion, which can account for some foot slippage. In the correction phase, forward-kinematic position measurements are used to correct the estimated state. This work was conducted on the StarlETH quadruped robot. If the stance feet orientations also remain constant, as is the case for many humanoids, Rotella et al. [68] showed that this ErEKF can be extended to allow forward-kinematic orientation measurements. The same group also formulated a similar unscented Kalman filter that uses forward-kinematic velocity measurements to correct inertial predictions and to accurately detect foot slip [16]. A detailed analysis of these filtering techniques combined with methods for incorporating computer vision can be found in Bloesch [17]. Due to the complexity involved in accurately formulating dynamic models, many groups have since adopted this IMU motion model approach to legged robot state estimation [84, 36, 46].

This combined inertial and kinematic filtering approach yields an estimate of the robot’s base pose and velocity. However, some legged robots require additional states to be estimated. Hwangbo et al. [50] formulated a probabilistic contact estimator for cases when contact sensors are unavailable. Xinjilefu et al. [84] developed a decoupled EKF that is able to estimate the full state of the humanoid robot ATLAS, including base states, joint angles, and joint velocities. Using proprioceptive sensing only, Bloesch et al. [15] proved that the absolute positions and yaw angles are unobservable. Thus, over time, estimates of these quantities will drift unboundedly. This is unacceptable for global mapping and planning algorithms; however, local elevation maps can still be obtained [37]. Fallon et al. [36] proposed a method for drift-free state estimation for the humanoid ATLAS. In their implementation, inertial and kinematic measurements were fused to yield accurate odometry. Point cloud data from a LiDAR sensor was used with a particle filter to localize the robot into a pre-built map. This approach provided corrections of position and yaw to obtain a drift-free estimate of the state. Nobili et al. [64] took a similar approach but used the Iterative Closest Point (ICP) algorithm to perform LiDAR-based point cloud matching. The algorithm was tested on the HyQ quadruped robot.

In this article, we propose using an InEKF to estimate the base pose and velocity of a general legged robot. The approach we take is most similar to Bloesch et al. [15], however we model the entire state as a single matrix Lie group as opposed to a decoupled state approach. This allows us to take advantage of the geometry and symmetry of the estimation problem to formulate autonomous error dynamics. In addition with our formulation, the linearizations are independent of the state estimate resulting in improved convergence properties, especially for poor state initializations. We formulate both world-centric and robo-centric state estimators highlighting the relation between the left- and right-invariant error dynamics. In addition, we provide exact analytical time-discretizations of both filters. The implemented filters can be run at high speeds ($>2000\mathop{\mathrm{Hz}}$) and can be directly used for accurate local odometry. We demonstrate this idea through a LiDAR terrain mapping application on a Cassie-series biped robot.

3 Theoretical background and preliminaries

We assume a matrix Lie group [43, 27] denoted $\mathcal{G}$ and its associated Lie Algebra denoted $\mathfrak{g}$. If elements of $\mathcal{G}$ are $n\times n$ matrices, then so are elements of $\mathfrak{g}$. When doing calculations, it is very convenient to let

[TABLE]

be the linear map that takes elements of the tangent space of $\mathcal{G}$ at the identity to the corresponding matrix representation so that the exponential map of the Lie group, $\exp:\mathbb{R}^{\mathrm{dim}\mathfrak{g}}\to\mathcal{G}$, is computed by

[TABLE]

where $\exp_{m}(\cdot)$ is the usual exponential of $n\times n$ matrices.

A process dynamics evolving on the Lie group with state at time $t$, $\textbf{X}_{t}\in\mathcal{G}$, is denoted by

[TABLE]

and $\bar{\textbf{X}}_{t}$ is used to denote an estimate of the state. The state estimation error is defined using right or left multiplication of $\textbf{X}_{t}^{-1}$ as follows.

Definition 1** (Left and Right Invariant Error).**

The right- and left-invariant errors between two trajectories $\textbf{X}_{t}$ and $\bar{\textbf{X}}_{t}$ are:

[TABLE]

where L is an arbitrary element of the group.

The following two theorems are the fundamental results for deriving an InEKF and show that by correct parametrization of the error variable, a wide range of nonlinear problems can lead to linear error equations.

Theorem 1** (Autonomous Error Dynamics [10]).**

A system is group affine if the dynamics, $f_{u_{t}}(\cdot)$, satisfies:

[TABLE]

for all $t>0$ and $\textnormal{{X}}_{1},\textnormal{{X}}_{2}\in\mathcal{G}$. Furthermore, if this condition is satisfied, the right- and left-invariant error dynamics are trajectory independent and satisfy:

[TABLE]

In the above, $\textbf{{I}}_{d}\in\mathcal{G}$ denotes the group identity element; to avoid confusion, we use I for a $3\times 3$ identity matrix, and $\textbf{I}_{n}$ for the $n\times n$ case. The following statements hold for both the left- and right-invariant error dynamics.

Define $\textbf{A}_{t}$ to be a $\mathrm{dim}\mathfrak{g}\times\mathrm{dim}\mathfrak{g}$ matrix satisfying

[TABLE]

For all $t\geq 0$, let $\boldsymbol{\xi}_{t}$ be the solution of the linear differential equation

[TABLE]

Theorem 2** (Log-Linear Property of the Error [10]).**

Consider the right-invariant error, $\boldsymbol{\eta}_{t}$, between two trajectories (possibly far apart). For arbitrary initial error $\boldsymbol{\xi}_{0}\in\mathbb{R}^{\mathrm{dim}\mathfrak{g}}$, if $\boldsymbol{\eta}_{0}=\exp(\boldsymbol{\xi}_{0})$, then for all $t\geq 0$,

[TABLE]

that is, the nonlinear estimation error $\boldsymbol{\eta}_{t}$ can be exactly recovered from the time-varying linear differential equation (3).

This theorem states that (3) is not the typical Jacobian linearization along a trajectory because the (left- or) right-invariant error on the Lie group can be exactly recovered from its solution. This result is of major importance for the propagation (prediction) step of the InEKF [10], where in the deterministic case, the log-linear error dynamics allows for exact covariance propagation.

Remark 1**.**

This indirect way of expressing the Jacobian of $g_{u_{t}}$ is from Barrau and Bonnabel [10]; it is used because we are working with a matrix Lie group viewed as an embedded submanifold of a set of $n\times n$ matrices.

During the correction step of a Kalman filter, the error is updated using incoming sensor measurements. If these observations take a particular form, then the linearized observation model and the innovation will also be autonomous [10]. This happens when the measurement, $\textbf{Y}_{t}$, can be written as either

[TABLE]

where b is a constant vector and $\textbf{V}_{t}$ is a vector of Gaussian noise.

The adjoint representation plays a key role in the theory of Lie groups and through this linear map we can capture the non-commutative structure of a Lie group.

Definition 2** (The Adjoint Map, see Hall [43] page 63).**

Let $\mathcal{G}$ be a matrix Lie group with Lie algebra $\mathfrak{g}$. For any $\textbf{X}\in\mathcal{G}$ the adjoint map, $\mathrm{Ad}_{\textbf{X}}:\mathfrak{g}\to\mathfrak{g}$, is a linear map defined as $\mathrm{Ad}_{\textbf{X}}({\boldsymbol{\xi}}^{\wedge})=\textbf{X}{\boldsymbol{\xi}}^{\wedge}\textbf{X}^{-1}$. Furthermore, we denote the matrix representation of the adjoint map by $\mathrm{Ad}_{\textbf{X}}$.

For more details on the material discussed above, along with the theory and proofs about the InEKF, we refer reader to Barrau [8], Barrau and Bonnabel [10, 11].

4 A motivating example: 3D orientation propagation

Suppose we are interested in estimating the 3D orientation of a rigid body given angular velocity measurements in the body frame, $\tilde{\boldsymbol{\omega}}_{t}\triangleq\mathrm{vec}(\omega_{x},\omega_{y},\omega_{z})\in\mathbb{R}^{3}$. This type of measurement can be easily obtained from a gyroscope.

There are several different parameterizations of $\mathrm{SO}(3)$; Euler angles, quaternions, and rotation matrices being the most common. If we let $\textbf{q}_{t}\triangleq\mathrm{vec}(q_{x},q_{y},q_{z})$ be a vector of Euler angles using the $R_{z}R_{y}R_{x}$ convention, then the orientation dynamics can be expressed as [33]

[TABLE]

Let $\delta\textbf{q}_{t}\triangleq\textbf{q}_{t}-\bar{\textbf{q}}_{t}\in\mathbb{R}^{3}$ be the error between the true and estimated Euler angles. The error dynamics can be written as a nonlinear function of the error variable, the inputs, and the state

[TABLE]

In order to propagate the covariance in an EKF, we need to linearize the error dynamics at the current state estimate, $\textbf{q}_{t}=\bar{\textbf{q}}_{t}$ (i.e. zero error). This leads to a linear error dynamics of the form:

[TABLE]

where $\bar{c}_{x}$, $\bar{s}_{x}$, and $\bar{t}_{x}$ are shorthand for $\cos(\bar{q}_{x})$, $\sin(\bar{q}_{x})$, and $\tan(\bar{q}_{x})$. The linear dynamics matrix clearly depends on the estimated angles, therefore bad estimates will affect the accuracy of the linearization and ultimately the performance and consistency of the filter.

Now instead, let’s use a rotation matrix to represent the 3D orientation, $\textbf{R}_{t}\in\mathrm{SO}(3)$. The dynamics can now be simply expressed as

[TABLE]

where $\left(\cdot\right)_{\times}$ denotes a $3\times 3$ skew-symmetric matrix. If we define the error between the true and estimated orientation as $\boldsymbol{\eta}_{t}\triangleq\textbf{R}_{t}^{\mathsf{T}}\bar{\textbf{R}}_{t}\in\mathrm{SO}(3)$, then the (left-invariant) error dynamics becomes

[TABLE]

Using this particular choice of state and error variable yields an autonomous error dynamics function (independent of the state directly). Since, $\mathrm{SO}(3)$ is a Lie Group, we can look at the dynamics of a redefined error that resides in the tangent space, $\boldsymbol{\eta}_{t}\triangleq\exp(\boldsymbol{\xi}_{t})$.

[TABLE]

After making a first-order approximation of the exponential map, the tangent space error dynamics become linear. In addition, this linear system only depends on the error, not the state estimate directly. In other words, wrong state estimates will not affect the accuracy of the linearization, which leads to better accuracy and consistency of the filter. For $\mathrm{SO}(3)$, this effect is well studied and has been leveraged to design the commonly used QEKFs333The set of quaternions, along with quaternion multiplication, actually forms a Lie group., [79, 75]. However, the extension to general matrix Lie groups, called the InEKF, was only recently developed by Barrau and Bonnabel [10].

In the above example, we utilized the first-order approximation for the exponential map of $\mathrm{SO}(3)$; $\exp(\boldsymbol{\xi}_{t})\approx\textbf{I}+\left(\boldsymbol{\xi}_{t}\right)_{\times}$. In general, one may ask how much accuracy is lost when making this approximation. The surprising result by Barrau and Bonnabel [10] is that this linearization is, in fact, exact. This is the basis of Theorem 2. If the initial error is known, the nonlinear error dynamics can be exactly recovered from this linear system.

Let’s demonstrate this theorem for the SO(3) example. Let $\boldsymbol{\eta}_{0}=\exp(\boldsymbol{\xi}_{0})$ be the initial left invariant error. We can show that $\boldsymbol{\eta}_{t}=\textbf{R}_{t}^{\mathsf{T}}\boldsymbol{\eta}_{0}\textbf{R}_{t}$ is a solution to the error dynamics equation (5) through differentiation.

[TABLE]

Once we replace the group error with the tangent space error and use the group’s adjoint definition to shift the rotation inside the exponential, we arrive at a simple expression for the tangent space error.

[TABLE]

Our final log-linear error dynamics can now be obtained by differentiating this new error expression.

[TABLE]

Again, this result shows that if the initial error is known, the nonlinear error dynamics can be exactly recovered from this linear system. In this work, we leverage these ideas to develop a contact-aided inertial observer for legged robots.

5 $\mathrm{SE}_{N+2}(3)$ Continuous Right-Invariant EKF

In this section, we derive a right-invariant extended Kalman filter (RIEKF) using IMU and contact motion models with corrections made through forward kinematic measurements. This RIEKF can be used to estimate the state of a robot that has an arbitrary (finite) number of points in contact with a static environment. While the filter is particularly useful for legged robots, the same theory can be applied for manipulators as long as the contact assumptions (presented in Section 5.2) are verified.

In order to be consistent with the standard InEKF theory, IMU biases are neglected for now. Section 7 provides a method for reintroducing the bias terms.

5.1 State Representation

As with typical aided inertial navigation, we wish to estimate the orientation, velocity, and position of the IMU (body) in the world frame [58, 38, 85]. These states are represented by $\textbf{R}_{\text{WB}}(t),{}_{\text{W}}\textbf{v}_{\text{B}}(t)$, and ${}_{\text{W}}\textbf{p}_{\text{WB}}(t)$ respectively. In addition, we append the position of all contact points (in the world frame), ${}_{\text{W}}\textbf{p}_{\text{WC}_{i}}(t)$, to the list of state variables. This is similar to the approach taken in Bloesch et al. [15], Bloesch [17].

The above collection of state variables forms a matrix Lie group, $\mathcal{G}$. Specifically, for $N$ contact points, $\textbf{X}_{t}\in\mathrm{SE}_{N+2}(3)$ can be represented by the following matrix:

[TABLE]

This Lie group is an extension of $\mathrm{SE}(3)$ and has been previously used to solve inertial navigation [8] and SLAM problems [19, 9]. In fact, the estimators derived in this work have a connection to the landmark-based SLAM problem. This connection is detailed later in Section 12.

Because the process and measurements models for each contact point, ${}_{\text{W}}\textbf{p}_{\text{WC}_{i}}(t)$, are identical, without loss of generality, we will derive all further equations assuming only a single contact point. Furthermore, for the sake of readability, we introduce the following shorthand notation:

[TABLE]

where the input $\textbf{u}_{t}$ is formed from the angular velocity and linear acceleration measurements coming from the IMU. It is important to note that these measurements are taken in the body (or IMU) frame. The Lie algebra of $\mathcal{G}$, denoted by $\mathfrak{g}$, is an $N+5$ dimensional square matrix. We use the hat operator, ${(\cdot)}^{\wedge}:\mathbb{R}^{3N+9}\to\mathfrak{g}$, to map a vector to the corresponding element of the Lie algebra. In the case of a single contact, for example, this function is defined by:

[TABLE]

where $\left(\cdot\right)_{\times}$ denotes a $3\times 3$ skew-symmetric matrix. The inverse operation is defined using the wedge operator, ${(\cdot)}^{\vee}:\mathfrak{g}\to\mathbb{R}^{3N+9}$. The matrix representation of the adjoint is given by:

[TABLE]

A closed form expression for the exponential map of $\mathrm{SE}_{N+2}(3)$ is given in Appendix B.

5.2 Continuous-Time System Dynamics

The IMU measurements are modeled as being corrupted by additive Gaussian white noise processes, per

[TABLE]

where $\mathcal{GP}$ denotes a Gaussian process and $\delta(t-t^{\prime})$ denotes the Dirac delta function. These are explicit measurements coming directly from a physical sensor. In contrast, the velocity of the contact point is implicitly inferred through a contact sensor; specifically, when a binary sensor indicates contact, the position of the contact point is assumed to remain fixed in the world frame, i.e. the measured velocity is zero. In order to accommodate potential slippage, the measured velocity is assumed to be the actual velocity plus white Gaussian noise, namely

[TABLE]

Using the IMU and contact measurements, the individual terms of the system dynamics can be written as:

[TABLE]

where g is the gravity vector and $\;\textit{{h}}_{R}(\tilde{\boldsymbol{\alpha}}_{t})$ is the measured orientation of the contact frame with respect to the IMU frame as computed through encoder measurements, $\tilde{\boldsymbol{\alpha}}_{t}\in\mathbb{R}^{M}$, and forward kinematics. Therefore, $\textbf{R}_{t}\;\textit{{h}}_{R}(\tilde{\boldsymbol{\alpha}}_{t})$ is a rotation matrix that transforms a vector from the contact frame to the world frame.

In matrix form, the dynamics can be expressed as

[TABLE]

with $\textbf{w}_{t}\triangleq\mathrm{vec}(\textbf{w}_{t}^{g},\;\textbf{w}_{t}^{a},\;\textbf{0}_{3,1},\;\textit{{h}}_{R}(\tilde{\boldsymbol{\alpha}}_{t})\textbf{w}_{t}^{v})$. The deterministic system dynamics, $f_{u_{t}}(\cdot)$, can be shown to satisfy the group affine property, (2). Therefore, following Theorem 1, the left- and right-invariant error dynamics will evolve independently of the system’s state.

Using Theorem 1, the right-invariant error dynamics is

[TABLE]

where the second term arises from the additive noise. The derivation follows the results in Barrau and Bonnabel [10] and is not repeated here.

Furthermore, Theorem 2 specifies that the invariant error satisfies a log-linear property. Namely, if $\textbf{A}_{t}$ is defined by $g_{u_{t}}(\exp(\boldsymbol{\xi}))\triangleq{(\textbf{A}_{t}\boldsymbol{\xi})}^{\wedge}+\mathcal{O}(||\boldsymbol{\xi}||^{2})$, then the log of the invariant error, $\boldsymbol{\xi}\in\mathbb{R}^{\mathrm{dim}\mathfrak{g}}$, approximately satisfies444With input noise, Theorem 2 no longer holds, and the linearization is only approximate. the linear system

[TABLE]

To compute the matrix $\textbf{A}_{t}$, we linearize the invariant error dynamics, $g_{u_{t}}(\cdot)$, using the first order approximation $\boldsymbol{\eta}_{t}^{r}=\exp(\boldsymbol{\xi}_{t})\approx\textbf{{I}}_{d}+{\boldsymbol{\xi}_{t}}^{\wedge}$ to yield

[TABLE]

With the above, we can express the prediction step of the RIEKF. The state estimate, $\bar{\textbf{X}}_{t}$, is propagated though the deterministic system dynamics, while the covariance matrix, $\textbf{P}_{t}$, is computed using the Riccati equation [61], namely,

[TABLE]

where the matrices $\textbf{A}_{t}$ and $\bar{\textbf{Q}}_{t}$ are obtained from (12) and (11),

[TABLE]

Remark 2**.**

For the right-invariant case, expression (14), $\textbf{A}_{t}$ is time-invariant and the time subscript could be dropped. However, in general it can be time-varying, therefore, we use $\textbf{A}_{t}$ throughout the paper.

5.3 Right-invariant Forward Kinematic Measurement Model

Let $\boldsymbol{\alpha}_{t}\in\mathbb{R}^{M}$ denote the vector of joint positions (prismatic or revolute) between the body and the contact point. We assume that the encoder measurements at time $t$ are corrupted by additive white Gaussian noise.

[TABLE]

Using forward kinematics, we measure the relative position of the contact point with respect to the body,

[TABLE]

where $\textbf{J}_{p}$ denotes the analytical Jacobian of the forward kinematics function. Using the state variables, the forward-kinematics position measurement becomes

[TABLE]

Re-written in matrix form, this measurement has the right-invariant observation form (4),

[TABLE]

Therefore, the innovation depends solely on the invariant error and the update equations take the form [10, Section 3.1.2]

[TABLE]

where $\exp(\cdot)$ is the exponential map corresponding to the state matrix Lie group, $\mathcal{G}$, $\textbf{L}_{t}$ is a gain matrix to be defined later, $\textbf{b}^{\mathsf{T}}=\begin{bmatrix}\textbf{0}_{1,3}&0&1&-1\end{bmatrix}$, and $\textbf{Y}_{t}^{\mathsf{T}}=\begin{bmatrix}\;\textit{{h}}_{p}^{\mathsf{T}}(\tilde{\boldsymbol{\alpha}}_{t})&0&1&-1\end{bmatrix}$. Because the last three rows of $\bar{\textbf{X}}_{t}\textbf{Y}_{t}-\textbf{b}$ are identically zero, we can express the update equations using a reduced dimensional gain, $\textbf{K}_{t}$, and an auxiliary selection matrix $\boldsymbol{\Pi}\triangleq\begin{bmatrix}\textbf{I}&\textbf{0}_{3,3}\end{bmatrix}$, so that $\textbf{L}_{t}\left(\bar{\textbf{X}}_{t}\textbf{Y}_{t}-\textbf{b}\right)=\textbf{K}_{t}\boldsymbol{\Pi}\bar{\textbf{X}}_{t}\textbf{Y}_{t}$ as detailed in Barrau [8].

Using the first order approximation of the exponential map, $\boldsymbol{\eta}_{t}^{r}=\exp(\boldsymbol{\xi}_{t})\approx\textbf{{I}}_{d}+{\boldsymbol{\xi}_{t}}^{\wedge}$, and dropping higher-order terms, we can linearize the update equation (18),

[TABLE]

Therefore,

[TABLE]

Taking ${(\cdot)}^{\vee}$ of both sides yields the linear update equation,

[TABLE]

Finally, we can write down the full state and covariance update equations of the RIEKF using the derived linear update equation and the theory of Kalman filtering [61, 2, 3] as

[TABLE]

where the gain $\textbf{K}_{t}$ is computed using

[TABLE]

and from (5.3), the matrices $\textbf{H}_{t}$ and $\bar{\textbf{N}}_{t}$ are given by

[TABLE]

5.4 Observability Analysis

Because the error dynamics are log-linear (c.f., Theorem 2), we can determine the unobservable states of the filter without having to perform a nonlinear observability analysis [8]. Noting that the linear error dynamics matrix in our case is time-invariant and nilpotent (with a degree of 3), the discrete-time state transition matrix is a polynomial in $\textbf{A}_{t}$,

[TABLE]

It follows that the discrete-time observability matrix is

[TABLE]

The last six columns (i.e., two matrix columns) of the observability matrix are clearly linearly dependent, which indicates the absolute position of the robot is unobservable. In addition, since the gravity vector only has a $z$ component, the third column of $\mathcal{O}$ is all zeros. Therefore, a rotation about the gravity vector (yaw) is also unobservable. This linear observability analysis agrees with the nonlinear observability results of Bloesch et al. [15], albeit with much less computation. Furthermore, as the error dynamics do not depend on the estimated state, there is no chance of the linearization spuriously increasing the numerical rank of the observability matrix [8]. This latter effect was previously known and studied by Bloesch et al. [15], and in order to resolve this problem, an observability-constrained EKF [49] was developed. In our proposed framework, by default, the discrete RIEKF has the same unobservable states as the underlying nonlinear system; hence, the developed discrete RIEKF intrinsically solves this problem.

6 Simulation Results

To investigate potential benefits or drawbacks of the proposed filter, we compare it against a state-of-the-art quaternion-based extended Kalman filter (QEKF), similar to those described by Bloesch et al. [15], Rotella et al. [68]. For implementation, the filter equations were discretized; see Appendix A for more details.

6.1 Quaternion-Based Filter Equations:

The choice of error variables is the main difference between the InEKF and the QEKF. Instead of the right-invariant error (1), a QEKF typically uses decoupled error states

[TABLE]

Using this definition of error, the QEKF deterministic error dynamics can be approximated as

[TABLE]

while the linearized observation matrix becomes

[TABLE]

The above linearizations are clearly dependent on the state estimate. Therefore, when the estimated state deviates from the true state, the linearizations are potentially wrong, reducing accuracy and consistency in the QEKF. In contrast, the deterministic right-invariant error dynamics are exactly log-linear (14). In addition, the linearized observation matrix for our InEKF (20) is also independent of the state estimate.

6.2 Convergence Comparison:

A dynamic simulation of a Cassie-series bipedal robot (described in Section 9) was performed in which the robot slowly walked forward after a small drop, accelerating from $0.0$ to $0.3\leavevmode\nobreak\ \mathop{\mathrm{m}}/\sec$. The discrete, simulated measurements were corrupted by additive white Gaussian noise, which are specified in Table 1 along with the initial state covariance values. The same values were used in both simulation and experimental convergence evaluations of the filters. The IMU bias estimation was turned off for these simulations. The simulation was performed with MATLAB and Simulink (Simscape Multibody™) where the simulation environment models ground contact forces with a linear force law (having a stiffness and damping term) and a Coulomb friction model. A typical walking gait is shown in Figure 2.

To compare the convergence properties of the two filters, 100 simulations of each filter were performed using identical measurements, noise statistics, initial covariance, and various random initial orientations and velocities. The initial Euler angle estimates were sampled uniformly from $-30\deg$ to $30\deg$. The initial velocity estimates were sampled uniformly from $-1.0\leavevmode\nobreak\ \mathop{\mathrm{m}}/\sec$ to $1.0\leavevmode\nobreak\ \mathop{\mathrm{m}}/\sec$. The pitch and roll estimates as well as the (body frame) velocity estimates for both filters are shown in Figure 3. Although both filters converge for this set of initial conditions, the proposed RIEKF converges considerably faster than the standard quaternion-based EKF.

6.3 Accuracy of Linearized Dynamics

The superior performance of the InEKF over the QEKF comes from the improved accuracy of the linearized error dynamics. As indicated by Theorem 2, the deterministic error dynamics of the InEKF are actually exact, while the QEKF version is only an approximation. To demonstrate this, a simulation was performed where propagation of the true error is compared to the propagation of the linearized error dynamics.

We first analyzed the deterministic dynamics. Given an initial error vector, $\boldsymbol{\xi}_{0}^{\text{true}}$, the initial state estimate for the InEKF was computed using the definition of right-invariant error (1), and the initial state estimate for the QEKF was computed using equation (21). The true state was initialized to the identity element. The true and estimated states for both filters were then propagated for 1 second (1000 time steps), using randomly sampled IMU measurements. The error states for both the InEKF and the QEKF were also propagated using their respective linearized error dynamics. The resulting error, $\boldsymbol{\xi}_{1}^{\text{true}}$, between the final estimated and true states were computed and compared to the propagated error states, $\boldsymbol{\xi}_{1}^{\text{prop}}$ to yield a measure of linearization accuracy, $\lVert\boldsymbol{\xi}_{1}^{\text{true}}-\boldsymbol{\xi}_{1}^{\text{prop}}\rVert$. This test was performed multiple times while linearly scaling the initial orientation error from $\mathrm{vec}(0,0,0)$ to $\mathrm{vec}(\pi/2,\pi/2,\pi/2)$. The results are shown in Figure 5.

As expected, when the initial error is zero, the difference between the true and propagated error is zero. This indicates that the linearized error dynamics for both the InEKF and the QEKF are correct. As the initial error increases, the difference between the true and propagated error states for the QEKF grows due to the decreased accuracy of the linearization. In contrast, the difference between the true and propagated error states for the InEKF are always exactly zero regardless of the initial error. In other words, assuming the initial error is known, the true propagated state can be exactly recovered from solving the linearized error dynamics system; see Theorem 2.

In the non-deterministic case (with sensor noise), Theorem 2 no longer holds. This can be seen in Figure 5, where the same test was performed, but with sensor noise corrupting the propagated state estimate. The difference between the true and propagated error is no longer exactly zero for the InEKF. However, the InEKF linearization remains more accurate due to the reduced sensitivity to initial state errors. This helps to further explain the improved convergence properties shown in Figure 3.

6.4 Covariance Ellipse Comparison

The error states in both the QEKF and the InEKF are assumed to be zero-mean Gaussian random vectors. However, due to the differing choice of error variables, the state uncertainty will differ. In the QEKF, all states and errors are decoupled (21). For example, the true position only depends on the position estimate and the position error, $\textbf{p}=\bar{\textbf{p}}+\delta\textbf{p}$. Therefore, the position estimate is a Gaussian centered at $\bar{\textbf{p}}$. In contrast, when using the InEKF, the position and orientation are actually coupled together, $\textbf{X}=\exp(\boldsymbol{\xi})\bar{\textbf{X}}$. Although $\boldsymbol{\xi}$ is a Gaussian random vector, after applying the group’s exponential map and matrix multiplication, the state estimate’s uncertainty distribution is no longer Gaussian. This distribution is known as a concentrated Gaussian on a Lie group [80, 81]. This type of distribution can often capture the underlying system uncertainty better than a standard Gaussian defined in Euclidean space [57, 7].

To demonstrate the difference, a simple simulation was performed where the Cassie robot walked forward for $8\leavevmode\nobreak\ \sec$ at an average speed of $1\leavevmode\nobreak\ m/s$. The standard deviation for the initial position uncertainty was set to $0.1\leavevmode\nobreak\ \mathop{\mathrm{m}}$ about each axis, while the standard deviation of the initial yaw uncertainty was set to $10\leavevmode\nobreak\ \deg$. A set of 10,000 particles sampled from this distribution were propagated forward to represent the robot’s true uncertainty distribution. After running both the InEKF and the QEKF, particles were sampled from the resulting filter covariances to provide a picture of the estimated position uncertainties. This result is shown in Figure 6.

The curved position distribution comes from the initial yaw uncertainty that continually grows due to its unobservability (Section 5.4). The InEKF is able to closely match this distribution since the samples are taken in the Lie algebra and are mapped to the group through the exponential map. This couples the orientation and position errors leading to a curved position distribution. In contrast, the QEKF position uncertainty can only have the shape of the standard Gaussian ellipse, which may not represent the true uncertainty well.

The InEKF can even accurately model the case of complete yaw uncertainty. To demonstrate this, the initial yaw standard deviation was set to $360\leavevmode\nobreak\ \deg$, and the same $8\leavevmode\nobreak\ \sec$ simulation was performed. Each ring in Figure 7 shows the sample position distribution spaced $2\leavevmode\nobreak\ \sec$ apart. This type of uncertainty cannot be captured with a standard Gaussian covariance ellipse. These examples illustrate that even if the means are identical, the covariance estimate of the InEKF can provide a more accurate representation of the state’s uncertainty than the standard QEKF.

In these examples, the initial covariance of all states (except yaw) was small. It is interesting to note what happens if we remove the dynamics noise and set the covariance for some states to exactly zero (rank-deficient covariance matrix). In this case, the initial covariance is supported by a subgroup and the InEKF will keep the state estimate within a time-dependent subset of the Lie group at all times. This theoretical result was proved by Chauchat et al. [25], Barrau and Bonnabel [12] to create EKFs with state equality constraints.

7 IMU bias augmentation

Implementation of an IMU-based state estimator on hardware typically requires modeling additional states, such as gyroscope and accelerometer biases. Unfortunately, as noted in Barrau [8], there is no Lie group that includes the bias terms while also having the dynamics satisfy the group affine property (2). Even though many of the theoretical properties of the RIEKF will no longer hold, it is possible to design an “imperfect InEKF” that still outperforms the standard EKF [8].

7.1 State Representation

The IMU biases are slowly varying signals that corrupt the measurements in an additive manner:

[TABLE]

These biases form a parameter vector that needs to be estimated as part of the RIEKF state,

[TABLE]

The model’s state now becomes a tuple of our original matrix Lie group and the parameter vector, $(\textbf{X}_{t},\boldsymbol{\theta}_{t})\in\mathcal{G}\times\mathbb{R}^{6}$. The augmented right-invariant error is now defined as

[TABLE]

Written explicitly, the right-invariant error is

[TABLE]

while the parameter vector error is given by

[TABLE]

7.2 System Dynamics

With IMU biases included, the system dynamics are now expressed as

[TABLE]

The IMU bias dynamics are modeled using the typical “Brownian motion” model, i.e., the derivatives are white Gaussian noise, to capture the slowly time-varying nature of these parameters,

[TABLE]

The deterministic system dynamics now depend on both the inputs, $\textbf{u}_{t}$, and the parameters, $\boldsymbol{\theta}_{t}$

[TABLE]

where $\bar{\boldsymbol{\omega}}_{t}\triangleq\tilde{\boldsymbol{\omega}}_{t}-\bar{\textbf{b}}_{t}^{g}$ and $\bar{\textbf{a}}_{t}\triangleq\tilde{\textbf{a}}_{t}-\bar{\textbf{b}}_{t}^{a}$ are the “bias-corrected” inputs. To compute the linearized error dynamics, the augmented right-invariant error (23) is first differentiated with respect to time,

[TABLE]

After carrying out the chain rule and making the first order approximation, $\boldsymbol{\eta}_{t}^{r}=\exp(\boldsymbol{\xi}_{t})\approx\textbf{{I}}_{d}+{\boldsymbol{\xi}_{t}}^{\wedge}$, the individual terms of the invariant error dynamics become

[TABLE]

Importantly, the augmented invariant error dynamics only depends on the estimated trajectory though the noise and bias errors, $\boldsymbol{\zeta}_{t}$ (this is expected because when there are no bias errors, there is no dependence on the estimated trajectory). A linear system can now be constructed from (27) to yield,

[TABLE]

where the noise vector is augmented to include the bias terms,

[TABLE]

7.3 Forward Kinematic Measurements

The forward kinematics position measurement (17) does not depend on the IMU biases. Therefore, the $\textbf{H}_{t}$ matrix can simply be appended with zeros to account for the augmented variables. The linear update equation becomes

[TABLE]

7.4 Final Continuous RIEKF Equations

The final “imperfect” RIEKF equations that include IMU biases can now be written down. The estimated state tuple is predicted using the following set of differential equations:

[TABLE]

The covariance of the augmented right invariant error dynamics is computed by solving the Riccati equation

[TABLE]

where the matrices $\textbf{A}_{t}$ and $\bar{\textbf{Q}}_{t}$ are now defined using (27),

[TABLE]

The estimated state tuple and its covariance are corrected though the update equations

[TABLE]

where the gains $\textbf{K}_{t}^{\boldsymbol{\xi}}$ and $\textbf{K}_{t}^{\boldsymbol{\zeta}}$ are computed from

[TABLE]

with the following measurement, output, and noise matrices,

[TABLE]

Remark 3**.**

The upper-right block of the new linearized dynamics matrix (28) is related to the adjoint of the current state estimate (7). Intuitively, this maps the bias error (measured in the body frame) to the world frame.

8 Addition and Removal of Contact Points

Sections 5 and 7 derived the equations for the RIEKF under the assumption that the contact point is unchanging with time. However, for legged robots, contacts are discrete events that are created and broken as a robot navigates through the environment. Therefore, it is important to be able to conveniently add and remove contact points states to and from the observer’s state.

8.1 Removing Contact Points

To remove a previous contact point from the state, we marginalize the corresponding state variable by simply removing the corresponding column and row from the matrix Lie group. The corresponding elements of the covariance matrix are also eliminated. This can be done through a simple linear transformation. For example, if the robot is going from one contact to zero contacts, then the newly reduced covariance would be computed by

[TABLE]

Remark 4**.**

This marginalization matrix, M, does not depend on the choice of right or left invariant error.

8.2 Adding Contact Points

When the robot makes a new contact with the environment, the state and covariance matrices need to be augmented. Special attention needs to be given to initialize the mean and covariance for the new estimated contact point. For example, if the robot is going from zero contacts to one contact, the initial mean is obtained though the forward kinematics relation

[TABLE]

In order to compute the new covariance, we need to look at the right-invariant error,

[TABLE]

Therefore, covariance augmentation can be done using the following linear map,

[TABLE]

Remark 5**.**

The error augmentation matrix, $\textbf{F}_{t}$, and the noise matrix, $\textbf{G}_{t}$, will depend on the choice of error variable. Here they are derived for the right invariant error case. The matrices will differ in the left invariant error formulation, as detailed in Section 10.

9 Experimental Results on Cassie Robot

We now present an experimental evaluation of the proposed contact-aided RIEKF observer using a 3D bipedal robot. The Cassie-series robot, shown in Figure 1, developed by Agility Robotics, has 20 degrees of freedom coming from the body pose, 10 actuators, and 4 springs. The robot is equipped with an IMU along with 14 joint encoders that can measure all actuator and spring angles. The proposed and baseline algorithms (along with the robot’s feedback controller) are implemented in MATLAB (Simulink Real-Time). The IMU (VectorNav-100) is located in the robot’s torso and provides angular velocity and linear acceleration measurements at $800\mathop{\mathrm{Hz}}$. The encoders provide joint angle measurements at $2000\mathop{\mathrm{Hz}}$. The robot has two springs on each leg that are compressed when the robot is standing on the ground. The spring deflections are measured by encoders and serve as a binary contact sensor. The controller used for these experiments was developed by Gong et al. [40].

9.1 Convergence Comparison

An experiment was performed where the robot walked forwards at approximately $0.3\leavevmode\nobreak\ \mathop{\mathrm{m}}/\sec$. The QEKF and the proposed RIEKF were run (off-line) 100 times using the same logged measurements, noise statistics, and initial covariance with random initial orientations and velocities. The noise statistics and initial covariance estimates are provided in Table 1. As with the simulation comparison presented in Section 6, the initial mean estimate for the Euler angles were uniformly sampled from $-30\deg$ to $30\deg$ and the initial mean estimate for velocities were sampled uniformly from $-1.0\leavevmode\nobreak\ \mathop{\mathrm{m}}/\sec$ to $1.0\leavevmode\nobreak\ \mathop{\mathrm{m}}/\sec$. Bias estimation was turned on and the initial bias estimate was obtained from processing the IMU data when the robot was static. The pitch and roll estimates as well as the (body frame) velocity estimates for both filters are shown in Figure 8. The “ground truth” trajectory estimates (black lines) were computed by initializing each filter with a good state estimate at a time before the beginning of the plot (to allow for convergence). The initial orientation was obtained from the VectorNav-100’s onboard EKF and the initial velocity was obtained through kinematics alone.

The experimental results for comparing filter convergence matches those of the simulation. The proposed RIEKF converges faster and more reliably in all 100 runs than the QEKF; therefore, due to the convenience of initialization and reliability for tracking the developed RIEKF is the preferred observer.

When the state estimate is initialized close to the true value, the RIEKF and QEKF have similar performance (black lines), because the linearization of the error dynamics accurately reflects the underlying nonlinear dynamics. However, when the state estimate is far from the true value, the simulation and experimental results show that RIEKF consistently converges faster than the QEKF. The relatively poor performance of the QEKF is due to the error dynamics being linearized around the wrong operating point, in which case the linear system does not accurately reflect the nonlinear dynamics. In addition, when bias estimation is turned off, the invariant error dynamics of the RIEKF do not depend on the current state estimate. As a result, the linear error dynamics can be accurately used even when the current state estimate is far from its true value, leading to better performance over the QEKF. Although this theoretical advantage is lost when bias estimation is turned on, the experimental results (shown in Figure 8) indicate that the RIEKF still is the preferred observer due to less sensitivity to initialization.

9.2 Motion Capture Experiment

In order to verify the accuracy of the InEKF state estimate, we performed a motion capture experiment in the University of Michigan’s M-Air facility. This outdoor space is equipped with 18 Qualisys cameras that allows for position tracking. We had the Cassie robot walk untethered for $60\leavevmode\nobreak\ \sec$ along an approximately 15 $\mathop{\mathrm{m}}$ path. A top-down view of the estimated trajectory is shown in Figure 9.

Although there is noticeable drift due to the unobservability of the position and yaw, the final position error accounts for less than $5\%$ of the distance traveled. This drift error is due to a combination of sensor noise and imperfect modeling of the robot’s kinematics which may introduce biases to the forward kinematic measurements. The orientation, velocity, and position estimates along with their $3\sigma$ covariance hulls are shown in Figure 10. Due to the unobservability of the yaw angle, the velocity estimate is given in the body frame instead of the world frame. The orientation is plotted using exponential coordinates, $\exp(\boldsymbol{\phi})=\textbf{R}$. Due to an inaccurate orientation estimate from the motion capture system, the “ground truth” for the orientation is given by the VectorNav-100, which runs a state-of-the-art QEKF that fuses angular velocity, linear acceleration, and magnetometer measurements to estimate orientation only.

As expected, the error for all observable states remains small. The absolute position and the orientation about the gravity vector are unobservable, so some drift will occur for these states. While the drift is largely imperceptible for $\boldsymbol{\phi}_{z}$, it is interesting to note that the covariance slowly grows over time due to this unobservability. More significant drift occurs on the absolute position states. This drift can be attributed to a combination of sensor noise, foot slip, and kinematic modeling errors.

In order to plot the $3\sigma$ covariance hull, the right-invariant error covariance needed to be converted to a covariance where the error is defined by Euclidean distance. Up to a first-order approximation, this mapping is done using:

[TABLE]

where the “Euclidean orientation error” is defined as $\delta\boldsymbol{\phi}_{t}\triangleq\boldsymbol{\phi}_{t}-\bar{\boldsymbol{\phi}}_{t}$, and the velocity and position errors match the QEKF error states (21). The matrix $\boldsymbol{\Gamma}_{1}(\bar{\boldsymbol{\phi}}_{t})$ is known as the left Jacobian of $\mathrm{SO}(3)$ and has an analytical form. Further explanation and the derivation of the above equation is given in Appendix C.

9.3 Long Odometry Experiment

In addition to providing accurate estimates of states vital for legged robot control (orientation and velocity), this InEKF can also provide reliable odometry for a higher-level mapping or SLAM system. To demonstrate the accuracy of long-term odometry, we had Cassie walk about $200\leavevmode\nobreak\ \mathop{\mathrm{m}}$ along a sidewalk around the University of Michigan’s Wave Field. In total, the walk took 7 minutes and 45 seconds. The estimated path from the InEKF overlaid onto Google Earth imagery is shown in Figure 11. A video of this experiment can be found at https://youtu.be/jRUltB_dMlo.

Even though the absolute position is unobservable, the odometry estimate from the InEKF contains low enough drift to keep the estimate on the sidewalk for the duration of the experiment. The final position estimate is within a few meters of the true position, and the yaw drift is imperceptible. This odometry estimate is readily available, as it only depends on inertial, contact, and kinematic data, which barring sensor failure, always exist. It does not require the use of any vision systems that may be susceptible to changes in environment or lighting conditions.

9.4 LiDAR Mapping Application

One application for the InEKF odometry is the building of local maps of the environment. We equipped the Cassie-series robot with a new torso that houses a Velodyne VLP-32C LiDAR. With the filter running, we can project each received packet of point cloud data into the world frame based on the current state estimate. This point cloud data can then be accumulated to create a map of the environment. Figure 12 shows a few still frames from several LiDAR mapping experiments. Videos of these results can be viewed at https://youtu.be/pNyXsZ5zVZk and https://youtu.be/nbQTQw0gJ-k.

10 Alternative Left-Invariant Formulation

For the derivations in Sections 5-8, we were assuming the use of the right-invariant error. This choice was due to the forward kinematic measurement having the right-invariant observation form. However, it is possible to derive a left-invariant form of this filter, which may be more appropriate to use when dealing with left-invariant observations. For example, GPS measurements are left-invariant observations for the world-centric observer; see Section 12. Written explicitly, the left-invariant error is

[TABLE]

After carrying out the chain rule and making the first order approximation, $\boldsymbol{\eta}_{t}^{l}=\exp(\boldsymbol{\xi}_{t})\approx\textbf{{I}}_{d}+{\boldsymbol{\xi}_{t}}^{\wedge}$, the individual terms of the left-invariant error dynamics become:

[TABLE]

Using these results, the log-linear left-invariant dynamics can be expressed using the following linear system

[TABLE]

where the dynamics and noise matrices are

[TABLE]

Similar to the right-invariant case, the dynamics only depend on the state through the IMU bias. When a left-invariant observation comes in, the state estimate is corrected using

[TABLE]

where the exponential map is now multiplied on the right side [10].

10.1 Switching Between Left and Right-Invariant Errors

Because forward kinematic measurements have the right invariant observation form, the innovation equations are only autonomous when the right invariant error is used. Fortunately, it is possible to switch between the left and right error forms through the use of the adjoint map.

[TABLE]

This transformation is exact, which means that we can easily switch between the covariance of the left and right invariant errors using

[TABLE]

Therefore, when handling a right-invariant observation, we can map the propagated left-invariant covariance to the right-invariant covariance temporarily, apply the right-invariant update equations (29), then map the corrected covariance back to the left-invariant form.

It is also possible to compute the log-linear left-invariant dynamics starting from the right-invariant form. Substituting (36) into the (right-invariant) covariance propagation equation (13) and solving for the left-invariant covariance yields

[TABLE]

Therefore, we learn that the right and left dynamics and noise matrices are related by the following expressions

[TABLE]

Remark 6**.**

Intuitively, the left-invariant error represents an error measured in the body frame of the robot, while the right-invariant error represents an error measured in the world or spatial frame. The frame of measurement dictates whether the exponential map appears on the right or left in the update equations. The error can be moved between these two frames using the adjoint map of the Lie group.

10.2 Adding New Contact Points

The process for removing a contact point from the state (marginalization) is identical to the right-invariant error case, described in Section 8.1. Likewise, when a new contact is detected, the state can be augmented using the same kinematics relation (31) as before. However, due to the change in error variable, the process for augmenting the covariance will be different.

In order to compute the new covariance, we need to look at the left-invariant error,

[TABLE]

Therefore, covariance augmentation can be done using the following linear map,

[TABLE]

11 Robo-centric Estimator

In this section, we derive a “robot-centric” version of the contact-aided InEKF where the estimated state is measured in the robot’s base (IMU) frame. When switching to a robot-centric model, the forward kinematics measurements take the left-invariant observation form. In addition, the left/right-invariant error dynamics equations are identical to the world-centric form, albeit swapped. The right-invariant error dynamics for the world-centric estimator are equivalent to the left-invariant error dynamics for the robo-centric estimator.

The properties of this filter are identical to the RIEKF derived in Section 5. However, in some cases this filter may be preferred as it directly estimates states that are useful for controlling a legged robot (namely the velocity measured in the body frame).

11.1 State and Dynamics

We are interested in estimating the same states as before, though measured in the robot’s body frame. Again, the state variables can form a matrix Lie group, $\mathcal{G}$. Specifically, for $N$ contact points, $\textbf{X}_{t}\in\mathrm{SE}_{N+2}(3)$ can be represented by the following matrix (which is simply the inverse of the world-centric state)555The negative sign on body velocity appears when inverting the world-centric state; $-{}_{\text{B}}\textbf{v}_{\text{B}}=-\textbf{R}_{\text{WB}}^{\mathsf{T}}\,{}_{\text{W}}\textbf{v}_{\text{B}}$. This sign is removed on the position vectors by swapping the start and end points, ${}_{\text{B}}\textbf{p}_{\text{BW}}=-\textbf{R}_{\text{WB}}^{\mathsf{T}}\,{}_{\text{B}}\textbf{p}_{\text{WB}}$.:

[TABLE]

Without loss of generality, assume a single contact point. Furthermore, for the sake of readability, we redefine our shorthand notation to be body-centric states,

[TABLE]

Using the bias corrected IMU measurements, the individual terms of the new robot-centric system dynamics can be derived as [16]

[TABLE]

Written in matrix form, this becomes

[TABLE]

with $\textbf{w}_{t}\triangleq\mathrm{vec}(\textbf{w}_{t}^{g},\,\textbf{w}_{t}^{a},\,\textbf{0}_{3,1},\;\textit{{h}}_{R}(\tilde{\boldsymbol{\alpha}}_{t})\textbf{w}_{t}^{v})$. The deterministic dynamics function $f_{u}(\cdot)$ can be shown to satisfy the group affine property (2). Therefore, the left (and right) invariant error dynamics depends solely on the invariant error.

We can derive the log-linear dynamics matrices for the body-centric estimator following a similar derivation process to the world-centric version in Section 5. In fact, without IMU bias, the linearization of the left-invariant dynamics for the body-centric estimator is the same as the right-invariant dynamics for the world-centric estimator;

[TABLE]

When IMU bias is included, the above relation still holds, though with the bias terms negated666If definition of bias error is negated, even these terms would remain the same.. The noise covariance matrices for the body-centric left/right-invariant propagation models are given by:

[TABLE]

which are also swapped versions of the world-centric noise matrices, after accounting for the redefinition of $\bar{\textbf{X}}_{t}$ as its inverse. A comparison of the world-centric and body-centric equations are given in Tables 2 and 3.

11.2 Left-Invariant Forward Kinematic Measurement Model

We use forward-kinematics to measure the relative position of the contact point with respect to the body, $\;\textit{{h}}_{p}(\tilde{\boldsymbol{\alpha}}_{t})$. Using the new robot-centric state variables, the measurement model (16) becomes

[TABLE]

Re-written in matrix form, this measurement will have now have the left-invariant observation form (4),

[TABLE]

The state update equation will take the left-invariant form (35), where the matrices $\textbf{H}_{t}$ and $\bar{\textbf{N}}_{t}$ can be derived to be:

[TABLE]

It is important to note that the forward kinematics measurement is a right-invariant observation for the world-centric estimator, while the same measurement becomes a left-invariant observation for the body-centric version. Also, the above linearization (46) is similar to the world-centric, right-invariant one (20). Namely, $\textbf{H}_{t}$ is simply negated, while $\bar{\textbf{N}}_{t}$ is identical after accounting for the redefinition of the state as its inverse.

12 Additional Sensor Observations

This article extends upon our conference paper results [46] where we originally presented the contact-aided InEKF. As described in earlier sections, this filter uses an inertial-contact dynamics model with corrections coming from forward kinematics. In addition to forward kinematics, other groups have discovered that a number of measurements common to robotics can also fit the invariant observation model (4). Bonnabel [18] developed an invariant observer that uses magnetometer and acceleration measurements to solve the attitude estimation problem. Barczyk and Lynch [4], Barrau [8] described methods for invariant observer design for GPS and magnetometer-aided navigation. Wu et al. [83] developed an InEKF to solve visual-inertial navigation. It has also been shown that an InEKF can be used for SLAM [19, 9, 86].

In particular, it is interesting to note the similarities between our contact-aided InEKF and landmark-based SLAM. In the simplest case, this SLAM problem involves jointly estimating the robots state along with the position of static landmarks in the environment. The robot is often assumed to have a sensor capable of measuring the position of the landmark relative to the robot. This formulation is identical to our developed InEKF with the contact positions acting as landmarks and forward kinematics measuring the relative translation between the base and contact frames. This similarity was also mentioned by Bloesch et al. [15]. However, there are a few notable differences. The contact frame velocity is assumed to be white noise to allow for foot slip, while landmarks are usually treated as static. Forward kinematics measurements often come at high frequencies ($2000\leavevmode\nobreak\ \mathop{\mathrm{Hz}}$ on Cassie). In contrast, landmarks measurements are often at a much lower frequency. Finally, with landmark observations, a data association problem often has to be solved which associates the measurement with a particular landmark state. This problem does not exist with forward kinematic measurements.

Due to the similarities between contacts and landmarks, in our state matrix, the contact states could be easily replaced with landmark states, ${}_{\text{W}}\textbf{p}_{\text{WL}_{i}}(t)$,

[TABLE]

or a combination of contact and landmark positions. In this way, it is possible to develop an observer that contains no unobservable states. Of course, like EKF-SLAM, the filter can become too computationally expensive to run in real-time if the number of landmarks grows too large.

Tables 2 and 3 give a summary of the left/right world-centric and robo-centric InEKF equations assuming a single contact and landmark position. The linearized observation matrix and observation type for several different sensors are also provided. In these tables, $\textbf{l}_{t}$ is shorthand for the true landmark position, and m denotes the true magnetic field vector. Using these tables, it is clear to see the relation between the left/right invariant error dynamics and the world/robo-centric formulations. All of these dynamics and observations are supported in an open-source C++ Library released alongside this article; available at: https://github.com/RossHartley/invariant-ekf.

13 Conclusion and future Work

Using recent results on a group-invariant form of the extended Kalman filter (EKF), this article derived an observer for a contact-aided inertial navigation system for a 3D legged robot. Contact and IMU sensors are available on all modern bipedal robots; therefore, the developed system has the potential to become an essential part of such platforms since an observer with a large basin of attraction can improve the reliability of perception and control algorithms. We also included IMU biases in the state estimator and showed that, while some of the theoretical guarantees are lost, in real experiments, the proposed system has better convergence performance than that of a commonly used quaternion-based EKF. Although the latter is a discrete EKF on a Lie group, it does not exploit symmetries present in the system dynamics and observation models, namely, invariance of the estimation error under a group action. In addition to the original right-invariant form of the filter, the left-invariant dynamics are discussed as well as a robo-centric version of the filter.

A series of experiments were conducted while running the filter on a Cassie-series biped robot. The accuracy of the invariant-EKF was demonstrated through a motion capture comparison, while long-term odometry was compared by overlaying the trajectory on Google Earth imagery. In addition to improving signals used for feedback control, the pose estimate from this filter can be used alongside a vision sensor to build maps of the environment. This was demonstrated on Cassie by building LiDAR-based point cloud maps while walking.

Future work includes developing an invariant smoother based on this filter. This would utilize the framework developed by Chauchat et al. [26] to potentially improve IMU preintegration [58, 38, 34] as well as contact preintegration [45] to perform SLAM. One interesting extension would be online estimation of kinematic parameters, which may help remove biases in the forward kinematic measurements. Another possibility for improvement of legged robot odometry is the detection of mode changes such as standing, flat-ground walking, and turning. During these modes, additional constraints may be inferred, which could improve state estimation [24]. Finally, additional experiments need to be conducted to explore the potential for visual-inertial-contact odometry using an InEKF as well as incorporating prior terrain information into the filter.

Appendix A Discretization of Filter Equations

In the preceding sections, the filters equations were in continuous time. However, in order to implement these filters using software and physical sensors, these equations need to be discretized. For our implementation, we assumed a zero-order hold on the inertial measurements, and performed analytical integration [34, 48]. In particular, analytical integration was important for the resulting error dynamics to satisfy Theorem 2.

A.1 World-centric Dynamics

This section demonstrates how to derive the deterministic, discrete time (world-centric) dynamics through analytical integration of the continuous state (24) and bias (25) dynamics. The contact, landmark, and bias dynamics are simply gaussian noise. Therefore, the discrete, deterministic dynamics are simply:

[TABLE]

Assuming a zero-order hold on the incoming IMU measurements between times $t_{k}$ and $t_{k+1}$, the orientation can be updated using the exponential map of $\mathrm{SO}(3)$:

[TABLE]

where $\Delta t\triangleq t_{k+1}-t_{k}$. Integrating the velocity dynamics yields an equation that involves the integral of the exponential map:

[TABLE]

Likewise, analytically solving for the discrete position dynamics involves computing the double integral:

[TABLE]

To solve these integrals, it is useful to define an auxiliary function [15]:

[TABLE]

which allows integrals to be easily expressed and computed using the taylor series form of the $\mathrm{SO}(3)$ exponential map.

[TABLE]

Closed form expressions also exist, allowing fast and easy computation of these quantities [75].

[TABLE]

Remark 7**.**

$\boldsymbol{\Gamma}_{0}(\boldsymbol{\phi})$ is simply the exponential map of $\mathrm{SO}(3)$, while $\boldsymbol{\Gamma}_{1}(\boldsymbol{\phi})$ is also known as the left Jacobian of $\mathrm{SO}(3)$ [27, 6].

Using these expressions, we can write down the discrete dynamics for the rotation, velocity, and positions states as:

[TABLE]

where the $t$ is dropped from the subscript for readability. These discrete dynamics are an exact integration of the continuous-time system under the assumption that the IMU measurements are constant over $\Delta t$.

A.2 Covariance Propagation

In order to propagate the covariance, a continuous-time Riccati equation needs to be solved.

[TABLE]

The analytical solution to the differential equation above is given by [61]:

[TABLE]

where the discrete noise covariance matrix is computed by

[TABLE]

and the state transition matrix, $\boldsymbol{\Phi}(t_{k+1},t_{k})$, satisfies

[TABLE]

The (world-centric) left-invariant error dynamics matrix only depends on the IMU inputs and the estimated bias terms, see Table 2. Since both are assumed to be constant between times $t_{k}$ and $t_{k+1}$, the state transition matrix can be simply computed from the matrix exponential.

[TABLE]

This state transition matrix also has an analytical solution of the form:

[TABLE]

where the individual terms are

[TABLE]

The matrices $\boldsymbol{\Psi_{1}}$ and $\boldsymbol{\Psi_{2}}$ involve computing the solution to a more complicated integral. However, these integrals still have analytical solutions which can be expressed easier after defining $\phi\triangleq\lVert\bar{\boldsymbol{\omega}}_{k}\rVert$ and $\theta\triangleq\phi\Delta t$.

[TABLE]

[TABLE]

The (world-centric) right-invariant error dynamics matrix depends on the the state estimates, $\bar{\textbf{R}}_{t}$, $\bar{\textbf{v}}_{t}$, and $\bar{\textbf{p}}_{t}$, which will change between times $t_{k}$ and $t_{k+1}$, see Table 2. Therefore, the state transition matrix will not simply be the matrix exponential, as in the left-invariant case. Solving (53) yields a state transition matrix of the form:

[TABLE]

where the individual terms can be analytically computed as

[TABLE]

Since the left/right-invariant errors are related through the adjoint, the two state transition and discrete noise matrices also satisfy the following relations [8].

[TABLE]

Therefore, the right-invariant state transition matrix can alternatively be computed using:

[TABLE]

which can simplify implementation since many software libraries already contain efficient methods to compute the matrix exponential.

Similar to the state transition matrices, the discrete noise covariance matrix (52) also has an analytical solution. In practice, this matrix is often approximated as:

[TABLE]

This approximated discrete noise matrix was used for all results in this article.

Appendix B Useful Lie Group Expressions

The matrix Lie group $\mathrm{SE}_{K}(3)$ is known as the group of $K$ direct isometries [8]. This group is comprised of a rotation matrix, $\textbf{R}\in\mathrm{SO}(3)$, and $K$ vectors in $\mathbb{R}^{3}$, $\textbf{p}_{1},\cdots,\textbf{p}_{K}$. Let X be an element of $\mathrm{SE}_{K}(3)$, which can be written as a $(3+K)\times(3+K)$ matrix:

[TABLE]

The group action is matrix multiplication. The adjoint is a linear map that can be used to move vectors between the tangent spaces of two group elements. The matrix representation of the adjoint of $\mathrm{SE}_{K}(3)$ is given by:

[TABLE]

The matrix representation of a vector, $\boldsymbol{\xi}\triangleq\mathrm{vec}(\boldsymbol{\phi},\boldsymbol{\xi}_{1},\cdots,\boldsymbol{\xi}_{K})\in\mathbb{R}^{3+3K}$, in the Lie algebra can obtained using the “hat” operator:

[TABLE]

where $\left(\boldsymbol{\phi}\right)_{\times}$ denotes the skew-symmetric matrix of a vector $\boldsymbol{\phi}=\mathrm{vec}(\boldsymbol{\phi}_{1},\boldsymbol{\phi}_{2},\boldsymbol{\phi}_{3})\in\mathbb{R}^{3}$ .

[TABLE]

The same vector can be moved to the Lie group through the exponential map:

[TABLE]

where $\boldsymbol{\Gamma}_{0}(\boldsymbol{\phi})$ is the exponential map of $\mathrm{SO}(3)$, and $\boldsymbol{\Gamma}_{1}(\boldsymbol{\phi})$ is the left Jacobian of $\mathrm{SO}(3)$.

Appendix C Error-state Conversions

Throughout this document, several versions of error states are used. These include the left/right invariant error (1), the QEKF error states (21), and the “Euclidean orientation error” used for plotting the covariance hull (33). It is often necessary to convert between these error states for plotting or when initializing the filters to provide fair comparisons. For example, if the QEKF and the InEKF are initialized with identical covariance matrices, the underlying distribution that they represent may be substantially different. This makes it impossible to show an accuracy comparison between the InEKF and QEKF with identical initial uncertainty. This section provides details on how to convert between these error states up to a first-order approximation.

When designing a QEKF, the orientation error can be defined in either the local or global frame [75]. These errors are equivalent to the left- and right-invariant errors for $\mathrm{SO}(3)$. In this document, the orientation error in the QEKF was chosen to be the error defined in the local frame (left-invariant error). Since the invariant errors are related through the group’s adjoint, the exact relation between right-invariant and QEKF orientation errors is

[TABLE]

When plotting the individual axes of the orientation error covariance hull, a “Euclidean orientation error” should be used. Let $\delta\boldsymbol{\phi}_{t}\triangleq\boldsymbol{\phi}_{t}-\bar{\boldsymbol{\phi}}_{t}$ be this Euclidean error where $\exp(\boldsymbol{\phi})\triangleq\textbf{R}$ is the exponential coordinate representation of a particular orientation. When the errors are small, a first-order approximation can be used to find a mapping between the right-invariant orientation error and this “Euclidean orientation error”.

[TABLE]

A similar first-order approximation can be used to find the relation between the right-invariant error and the QEKF velocity errors.

[TABLE]

The same process can be repeated for the position states.

[TABLE]

Acknowledgments

The authors would like to thank Yukai Gong for the development of the feedback controller utilized in the Cassie experiments as well as Bruce Huang, Ray Zhang, Lu Gan, Zhenyu Gan, Omar Harib, Eva Mungai, and Grant Gibson for their help in collecting experimental data. Funding for R. Hartley and M. Ghaffari is given by the Toyota Research Institute (TRI), partly under award number N021515, however this article solely reflects the opinions and conclusions of its authors and not TRI or any other Toyota entity. Funding for J. Grizzle was in part provided by TRI and in part by NSF Award No. 1808051.

IMU inertial measurement unit EKF extended Kalman filter ErEKF error-state extended Kalman filter QEKF quaternion-based extended Kalman filter InEKF invariant extended Kalman filter RIEKF right-invariant extended Kalman filter LIEKF lef- invariant extended Kalman filter INS inertial navigation system FK forward kinematic DOF degrees of freedom SLAM simultaneous localization and mapping ICP Iterative Closest Point