Learning-Assisted Multi-IMU Proprioceptive State Estimation for Quadruped Robots

1. Introduction

Accurate state estimation is critical for balance and motion planning of quadruped robots on challenging terrains. The requirements of real-time implementation and high accuracy for state estimation in the real world restrict the sensor and algorithm selection [1].

Exteroception and proprioception are two distinct sensory modalities that play crucial roles in state estimation for legged robots [2]. While exteroception has been investigated for state estimation with external environment information from sensors such as cameras [3] and lidars [4], exteroceptive sensors may be prone to environmental disturbances and perceptual limitations (e.g., occlusions, varying lighting conditions), and provide large amounts of high-dimension sensory data that can be computationally demanding for real-time processing, which limits its applicability. Instead, proprioception involves sensing the internal state of the robot, including internal sensors like joint encoders, force sensors, and inertial measurement units (IMUs) that provide immediate feedback about the robot’s position, orientation, and internal forces in real time. Proprioceptive state estimation has become a popular design due to its low latency, high energy efficiency, and low computational demands [5]. The challenges of proprioceptive state estimation lie in calibration, drift, and sensor noise.

Extended Kalman filter (KF) is a commonly used approach for sensor fusion and state estimation [6,7]. EKF integrates measurements from multiple sensors (e.g., encoders, IMUs) with the robot’s motion model to estimate the state variables while considering measurement noise and uncertainties. Bloesch et al. [8] presented a consistent fusion of leg kinematics and IMU using an observability-constrained EKF for state estimation. Furthermore, Yang et al. [9] proposed online kinematic calibration for legged robots to reduce the velocity estimation errors caused by inaccurate leg-length knowledge. Then, Multi-IMU Proprioceptive Odometry (MIPO) [10] was developed to improve upon KF-based Proprioceptive Odometry (PO) methods that only use a single IMU. It is noted that MIPO focuses on PO rather than the full state estimation and uses the mocap system for measuring the body yaw. However, EKF is capable of diverging in some challenging situations. Invariant EKF (IEKF) has been increasingly investigated to enhance convergence of state estimation. The state is defined on a Lie group and dynamics satisfy a particular group-affine property in IEKF. The invariance of the estimation error with respect to a Lie group action leads to the estimation error satisfying a log-linear autonomous differential equation on the Lie algebra, allowing the design of a nonlinear state estimator with strong convergence properties. Hartley et al. [11] derived an IEKF for a system containing IMU and contact sensor dynamics, with forward kinematic correction measurements. However, few works consider the IEKF for multi-IMU proprioceptive state estimation in the context of quadruped robots. Xavier et al. [12] considered multiple IMUs distributed throughout the robot’s structure for state estimation of humanoid robots. In this work, we extend IEKF for state estimation of a quadruped robot with a body IMU and four foot IMUs, which makes a good benchmark.

Moreover, machine learning has been employed to address the challenges of proprioceptive state estimation. In particular, learning-based approaches have been developed for reliable contact detection and measurements estimation with reduced sensor demands under challenging scenarios. Lin et al. [13] developed a learning-based contact estimator for IEKF to bypass the need for physical contact sensors. Buchanan et al. [14] proposed a learned inertial displacement measurement to improve state estimation in challenging scenarios where leg odometry is unreliable, such as slipping and compressible terrains. Teng et al. [15] proposed a neural measurement network (NMN) to estimate the contact probabilities of feet and the body linear velocity for IEKF. The inputs to the NMN consist of current body acceleration, body angular velocity, joint angles, joint velocities, and previous positional joint targets. Moreover, Liu et al. [16] employed a physics-informed neural network in conjunction with an unscented KF using proprioceptive sensory data to enhance the state estimation process. Inspired by the previous works, this paper aims to develop an accurate full-state estimator using only proprioceptive information. In particular, built upon the MIPO, we use a 1D CNN learned from multi-IMU and other internal sensor data with groundtruth body orientations to estimate the body yaw, to avoid the mocap or sophisticated visual-inertial-leg odometry algorithm (e.g., [17]) for accurate yaw measurements/estimates. Then, the yaw angles predicted by the CNN model are used by EKF to estimate states. Furthermore, we compare this learning-assisted EKF approach with the IEKF using multi-IMU measurements.

The main contribution of this paper lies in presenting a learning-assisted multi-IMU proprioceptive state estimation (LMIPSE) approach that provides accurate full-state estimation using only proprioception. The remainder of the paper is organized as follows: Section 2 provides the problem formulation and related preliminaries; Section 3 presents the right-invariant EKF (RI-EKF) using multi-IMU measurements; Section 4 presents the proposed LMIPSE; Section 5 provides experimental results; and finally, Section 6 summarizes this paper.

2. Problem Formulation and Related Preliminaries

The states to estimate for a quadruped robot include the orientation $\mathbf{R}\in SO\left(3\right)$, velocity $\mathbf{v}\in {\mathbb{R}}^3$, and position $\mathbf{p}\in {\mathbb{R}}^3$ of the body in the world frame W. The sensor measurements considered for updating state estimations in this paper are from 5 IMUs and 4 joint encoders. We use $I_0$ to denote the body IMU frame which is aligned with the body frame B and use $\left\{I_i\right|i=1,2,3,4\}$ to denote the frames of the IMUs attached to the legs. Additionally, $C_i$ denotes the frames of the foot-end contact points and $i=1,2,3,4$ correspond to the front left, front right, hind left, and hind right foot, respectively.

2.1. Measurements and System Model

Consider the IMU measurements model described by

$$\begin{array}{cc}\hfill {\tilde{\omega }}_{i,t}& ={\omega }_{i,t}+{\mathbf{b}}_{i,t}^g+{\mathbf{w}}_{i,t}^g,\hfill \end{array}$$

(1a)

$$\begin{array}{cc}\hfill {\tilde{\mathbf{a}}}_{i,t}& ={\mathbf{a}}_{i,t}+{\mathbf{b}}_{i,t}^a+{\mathbf{w}}_{i,t}^a,\hfill \end{array}$$

(1b)

$$\begin{array}{cc}\hfill {\dot{\mathbf{b}}}_{i,t}^g& ={\mathbf{w}}_{i,t}^{b,g},{\dot{\mathbf{b}}}_{i,t}^a={\mathbf{w}}_{i,t}^{b,a},\hfill \end{array}$$

(1c)

where ${\tilde{\omega }}_{i,t}$, ${\mathbf{b}}_{i,t}^g$, and ${\mathbf{w}}_{i,t}^g$ denote the measurement, bias, and white Gaussian noise of the gyroscope in the i-th IMU, respectively; ${\tilde{\mathbf{a}}}_{i,t}$, ${\mathbf{b}}_{i,t}^a$, and ${\mathbf{w}}_{i,t}^a$ denote the measurement, bias, and white Gaussian noise of the accelerometer in the i-th IMU, respectively; Eq. (1c) considers random walk model of the bias terms of the IMU, and ${\mathbf{w}}_{i,t}^{b,g}$ and ${\mathbf{w}}_{i,t}^{b,a}$ are white Gaussian noise. Given the angular velocity ${\tilde{\omega }}_{0,t}$ and linear acceleration ${\tilde{\mathbf{a}}}_{0,t}$ measured by the body IMU as well as the contacts to the ground and the joint positions ${\tilde{\alpha }}_{i,t}$ from the joint encoders, the system model can be described by

$$\begin{array}{cc}\hfill {\dot{\mathbf{R}}}_t& ={\mathbf{R}}_t{({\tilde{\omega }}_{0,t}-{\mathbf{b}}_{0,t}^g-{\mathbf{w}}_{0,t}^g)}_{\times },\hfill \end{array}$$

(2a)

$$\begin{array}{cc}\hfill {\dot{\mathbf{v}}}_t& ={\mathbf{R}}_t({\tilde{\mathbf{a}}}_{0,t}-{\mathbf{b}}_{0,t}^a-{\mathbf{w}}_{0,t}^a)+\mathbf{g},\hfill \end{array}$$

(2b)

$$\begin{array}{cc}\hfill {\dot{\mathbf{p}}}_t& ={\mathbf{v}}_t,\hfill \end{array}$$

(2c)

$$\begin{array}{cc}\hfill {\dot{\mathbf{d}}}_{i,t}& ={\mathbf{R}}_t{\mathbf{h}}_R\left({\tilde{\alpha }}_{i,t}\right)(\mathbf{0}-{\mathbf{w}}_{i,t}^v),\hfill \end{array}$$

(2d)

where ${(·)}_{\times }$ denotes a skew symmetric matrix; $\mathbf{g}$ is the gravity vector; ${\mathbf{d}}_{i,t}$ is the position of the i-th contact point in the world frame; ${\mathbf{h}}_R\left({\alpha }_t\right)=-{\mathbf{J}}_i\dot{\alpha }-{\left(\omega \right)}_{\times }{\mathbf{r}}_i\left(\alpha \right)$ with ${\mathbf{J}}_i={\displaystyle \frac{\partial {\mathbf{r}}_i}{\partial \alpha }}$ and ${\mathbf{r}}_i$ denoting the location of the i-th contact point in the body frame. It is noted that Eq. (2d) assumes zero measured velocity for the contact point.

2.2. Right-Invariant EKF

The right-invariant EKF (RI-EKF) has proven consistent in orientation estimation by maintaining the structure of the orientation space and effective in drift mitigation by properly propagating uncertainties through the state space [18]. The RI-EKF for one body IMU and the contact process model with a forward kinematics measurement model has been derived in [11,19]. The prediction step and the update step of the RI-EKF are as follows.

2.2.1. Prediction Step

Assuming zero-order hold to the input and performing Euler integration from time $t_k$ to $t_{k+1}$, the discrete dynamics are described by

$$\begin{array}{cc}\hfill {\widehat{\mathbf{R}}}_{k+1}^{-}& ={\widehat{\mathbf{R}}}_k^{+}exp(({\tilde{\omega }}_k-{\left({\widehat{\mathbf{b}}}_{0,k}^g\right)}^{+})\Delta t),\hfill \end{array}$$

(3a)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{v}}}_{k+1}^{-}& ={\widehat{\mathbf{v}}}_k^{+}+{\widehat{\mathbf{R}}}_k^{+}({\tilde{\mathbf{a}}}_k-{\left({\widehat{\mathbf{b}}}_{0,k}^a\right)}^{+})\Delta t+\mathbf{g}\Delta t,\hfill \end{array}$$

(3b)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{p}}}_{k+1}^{-}& ={\widehat{\mathbf{p}}}_k^{+}+{\widehat{\mathbf{v}}}_k^{+}\Delta t+{\displaystyle \frac{1}{2}}{\widehat{\mathbf{R}}}_k^{+}({\tilde{\mathbf{a}}}_k-{\left({\widehat{\mathbf{b}}}_{0,k}^a\right)}^{+})\Delta t^2+{\displaystyle \frac{1}{2}}\mathbf{g}\Delta t^2,\hfill \end{array}$$

(3c)

$$\begin{array}{cc}\hfill {\widehat{\mathbf{d}}}_{k+1}^{-}& ={\widehat{\mathbf{d}}}_k^{+},{\left({\widehat{\mathbf{b}}}_{0,k+1}^g\right)}^{-}={\left({\widehat{\mathbf{b}}}_{0,k}^g\right)}^{+},{\left({\widehat{\mathbf{b}}}_{0,k+1}^a\right)}^{-}={\left({\widehat{\mathbf{b}}}_{0,k}^a\right)}^{+},\hfill \end{array}$$

(3d)

$$\begin{array}{cc}\hfill P_{k+1}& ={\Phi }_k{P}_k{\Phi }_k^{\top }+Q_k,\hfill \end{array}$$

(3e)

where $\Delta t=t_{k+1}-t_k$, ${(·)}_k^{+}$ denotes the corrected estimation with all the measurements until $t_k$ and ${(·)}_{k+1}^{-}$ denotes the one-step-ahead prediction using the dynamic model with the corrected estimates at $t_k$; ${\Phi }_k={exp}_{\mathrm{m}}(A_k\Delta t)$ is the discrete state transformation matrix with $A_k{:=A|}_{t=t_k}$ and

$$A=\left[\begin{array}{cccccc}{\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -\widehat{\mathbf{R}}& {\mathbf{0}}_3\\ {\mathbf{g}}_{\times }& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\left(\widehat{\mathbf{v}}\right)}_{\times }\widehat{\mathbf{R}}& -\widehat{\mathbf{R}}\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\left(\widehat{\mathbf{p}}\right)}_{\times }\widehat{\mathbf{R}}& {\mathbf{0}}_3\\ {\mathbf{0}}_3& I_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\left(\widehat{\mathbf{d}}\right)}_{\times }\widehat{\mathbf{R}}& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3\end{array}\right];$$

being the system matrix of the linearized right-invariant error dynamics; the estimation error $\xi ={[{\xi }_R^{\top },{\xi }_v^{\top },{\xi }_p^{\top },{\xi }_d^{\top },e_{b_0^g}^{\top },e_{a_0^g}^{\top }]}^{\top }$ with $e_{b_0^g}={\widehat{\mathbf{b}}}_0^g-{\mathbf{b}}_0^g$ and $e_{b_0^a}={\widehat{\mathbf{b}}}_0^a-{\mathbf{b}}_0^a$; $Q_k\approx {\Phi }_k\widehat{Q}{\Phi }_k^{\top }\Delta t$ with $\widehat{Q}=\left[\begin{array}{cc}{\mathrm{Ad}}_{{\widehat{X}}_t}& {\mathbf{0}}_{12\times 6}\\ {\mathbf{0}}_{6\times 12}& I_6\end{array}\right]\mathrm{Cov}\left({\mathbf{w}}_t\right){\left[\begin{array}{cc}{\mathrm{Ad}}_{{\widehat{X}}_t}& {\mathbf{0}}_{12\times 6}\\ {\mathbf{0}}_{6\times 12}& I_6\end{array}\right]}^{\top }$ and

$${\mathbf{w}}_t:=\mathrm{vec}({\mathbf{w}}_t^g,{\mathbf{w}}_t^a,{\mathbf{0}}_{3\times 1},{\mathbf{h}}_R\left({\tilde{\alpha }}_t\right){\mathbf{w}}_t^v,{\mathbf{w}}_{0,t}^{b,g},{\mathbf{w}}_{0,t}^{b,a});$$

${exp}_{\mathrm{m}}(·)$ is the matrix exponential; ${\mathrm{Ad}}_X\left({\xi }^{\wedge }\right)=X{\xi }^{\wedge }{X}^{-1}$ denotes the adjoint map and ${\xi }^{\wedge }$ is the associated element in a $n\times n$ matrix Lie group for $\xi \in {\mathbb{R}}^n$.

2.2.2. Update Step

$$\begin{array}{cc}\hfill {\widehat{X}}^{+}& =exp\left({\delta }_{IMU}\right){\widehat{X}}^{-},\hfill \end{array}$$

(4a)

$$\begin{array}{cc}\hfill {\mathbf{b}}_g^{+}& ={\mathbf{b}}_g^{-}+{\delta }_{bg},{\mathbf{b}}_a^{+}={\mathbf{b}}_a^{-}+{\delta }_{ba},\hfill \end{array}$$

(4b)

$$\begin{array}{cc}\hfill P^{+}& =(I-KH)P^{-}{(I-KH)}^{\top }+KN{K}^{\top },\hfill \end{array}$$

(4c)

where

$$X_t:=\left[\begin{array}{cccc}{\mathbf{R}}_t& {\mathbf{v}}_t& {\mathbf{p}}_t& {\mathbf{d}}_t\\ {\mathbf{0}}_{1\times 3}& 1& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 1& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 1\end{array}\right]$$

(5)

is the embedding of the states in the matrix Lie group1; $\delta :={[{\delta }_{IMU}^{\top },{\delta }_{bg}^{\top },{\delta }_{ba}^{\top }]}^{\top }=K\Pi \left(\widehat{X}Y\right)$ is the correction term with K being the Kalman gain and $\Pi =[I_3,{\mathbf{0}}_3]$ for selecting the first three rows of the right-hand-side matrix. The Kalman gain is given by

$$\begin{array}{cc}\hfill N& =\widehat{\mathbf{R}}{\mathbf{J}}_v\left(\tilde{\alpha }\right)\mathrm{Cov}\left({\mathbf{w}}_{\alpha }\right){\mathbf{J}}_v^{\top }\left(\tilde{\alpha }\right){\widehat{\mathbf{R}}}^{\top }+Q_{\alpha },\hfill \end{array}$$

(6a)

$$\begin{array}{cc}\hfill H& =[{\mathbf{0}}_3,{\mathbf{0}}_3,-I_3,I_3,{\mathbf{0}}_3,{\mathbf{0}}_3],\hfill \end{array}$$

(6b)

$$\begin{array}{cc}\hfill S& =H{P}^{-}{H}^{\top }+N,\hfill \end{array}$$

(6c)

$$\begin{array}{cc}\hfill K& =P{H}^{\top }{S}^{-1},\hfill \end{array}$$

(6d)

where $Q_{\alpha }$ is the covariance matrix of the additive white Gaussian noise from the leg kinematics measurements.

2.3. Multi-IMU Proprioceptive Odometry

Multiple IMUs, including body and foot IMUs, have been used with joint encoders for legged robots to achieve low-drift long-term position and velocity estimation using an EKF. The foot IMUs can address the limitations of the PO using the single body IMU by enabling the updates of the foot positions during the airborne phase and considering the rolling during the contact without the zero-velocity assumption for the Leg Odometry velocity. In particular, foot IMU data was employed to update foot velocities in the prediction step and determine foot contact models and slip for the measurement model [10]. The process model for MIPO is

$${\widehat{x}}_{k+1}=\left[\begin{array}{c}{\widehat{\theta }}_{k+1}\\ {\widehat{\mathbf{v}}}_{k+1}\\ {\widehat{\mathbf{p}}}_{k+1}\\ {\widehat{\mathbf{d}}}_{i,k+1}\\ {\dot{\widehat{\mathbf{d}}}}_{i,k+1}\end{array}\right]=\left[\begin{array}{c}{\widehat{\theta }}_k+(\Omega \left({\widehat{\theta }}_k\right){(}^B\omega ))\Delta t\\ {\widehat{\mathbf{v}}}_k+(\mathbf{R}\left({\widehat{\theta }}_k\right){\mathbf{a}}_{\mathbf{0}}-\mathbf{g})\Delta t\\ {\widehat{\mathbf{p}}}_k+{\widehat{\mathbf{v}}}_k\Delta t\\ {\widehat{\mathbf{d}}}_{i,k}+{\dot{\widehat{\mathbf{d}}}}_{i,k}\Delta t\\ {\dot{\widehat{\mathbf{d}}}}_{i,k}+(\mathbf{R}\left({\widehat{\theta }}_k\right){\mathbf{h}}_R\left(\alpha \right){\mathbf{a}}_{\mathbf{i}}-\mathbf{g})\Delta t\end{array}\right]$$

(7)

where $\Omega \left({\widehat{\theta }}_k\right)$ originates from the derivative of $\theta ={\left[\begin{array}{ccc}{\theta }_r& {\theta }_p& {\theta }_y\end{array}\right]}^{\top }$ with ${\theta }_r$, ${\theta }_p$, $\theta y$ denoting the roll, pitch, and yaw angles, respectively. The measurement model is

$$h({\widehat{x}}_k,\alpha ,{\omega }_i)=\left[\begin{array}{c}\mathbf{R}{\left({\widehat{\theta }}_k\right)}^{\top }(\widehat{\mathbf{d}}-{\widehat{\mathbf{p}}}_k)\\ \mathbf{R}{\left({\widehat{\theta }}_k\right)}^{\top }({\widehat{\mathbf{v}}}_{i,k}-{\widehat{\mathbf{v}}}_k)\\ {\dot{\widehat{\mathbf{d}}}}_{i,k}{-}^W{\omega }_i\times \mathbf{s}\end{array}\right]$$

(8)

where the last term is based on the pivoting model [10]; × denotes the cross product of two vectors; ${}^W{\omega }_i={\mathbf{h}}_R\left(\tilde{\alpha }\right){\omega }_i$; using $d_0$ to denote the distance between the foot center and the foot surface and $\mathbf{n}=\mathbf{R}\left({\widehat{\theta }}_k\right){\mathbf{h}}_R\left(\tilde{\alpha }\right)$ to denote the contact normal vector expressed in the world frame, $\mathbf{s}(\widehat{x},\tilde{\alpha })=-s_0·\mathbf{n}/\parallel \mathbf{n}\parallel$ is the pivoting vector pointing from the contact point to the body center. The last term in (8) can be used for the update step only when a foot is in contact with the ground. To detect the contact and slip, a statistical test based on Mahalanobis distance is used in [10] as follows:

$$\sqrt{z^T{S}_z^{-1}z}<\sigma ,$$

(9)

where $z={\dot{\widehat{\mathbf{d}}}}_{i,k}{-}^W{\omega }_i\times \mathbf{s}$ and $S_z$ is the covariance matrix of z estimated by the EKF; $\sigma$ is the hypothesis testing threshold. The foot is recognized as being in non-slipping contact if (9) is satisfied. Since RI-EKF shows superior performance, we will present the RI-EKF approach for state estimation using multi-IMU measurements.

3. RI-EKF Using Multi-IMU Measurements

In this section, we present RI-EKF using multi-IMU measurements as a strong benchmark. Considering the linear velocity of the foot as a state in the system dynamics, the augmented robot state $\overline{x}=[\theta ,{\mathbf{v}}_0,\mathbf{p},{\mathbf{d}}_i,{\mathbf{v}}_i]$ where ${\mathbf{v}}_i$ denotes the linear velocity of the foot i contact foot in the world frame and ${\dot{\mathbf{v}}}_{i,k}=\mathbf{R}\left({\theta }_k\right){\mathbf{h}}_R\left(\alpha \right){\mathbf{a}}_i+\mathbf{g}$. Then, we consider

$${\overline{X}}_t:=\left[\begin{array}{ccccc}{\mathbf{R}}_t& {\mathbf{v}}_{0,t}& {\mathbf{p}}_t& {\mathbf{d}}_{i,t}& {\mathbf{v}}_{i,t}\\ {\mathbf{0}}_{1\times 3}& 1& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 1& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 1& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 1\end{array}\right].$$

The foot IMUs can provide the estimation of the foot velocities and accelerations even during slipping by

$$\begin{array}{cc}\hfill {\dot{d}}_{i,t}& =v_{i,t},\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill {\dot{v}}_{i,t}& ={\mathbf{R}}_t{\mathbf{h}}_R\left(\tilde{\alpha }\right)({\tilde{\mathbf{a}}}_{i,t}-{\mathbf{b}}_{i,t}^a-{\mathbf{w}}_{i,t}^a)+\mathbf{g}.\hfill \end{array}$$

(11)

Consequently, we have

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\overline{X}}_t& =\left[\begin{array}{ccccc}{\mathbf{R}}_t{\left({\tilde{\omega }}_{0,t}\right)}_{\times }& {\mathbf{R}}_t{\tilde{\mathbf{a}}}_{0,t}+\mathbf{g}& {\mathbf{v}}_{0,t}& {\mathbf{0}}_{3\times 1}& {\mathbf{R}}_t{\mathbf{h}}_R\left(\tilde{\alpha }\right){\tilde{\mathbf{a}}}_{i,t}+\mathbf{g}\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\end{array}\right]\hfill \\ & -\left[\begin{array}{ccccc}{\mathbf{R}}_t& {\mathbf{v}}_{0,t}& {\mathbf{p}}_t& {\mathbf{d}}_{i,t}& {\mathbf{v}}_{i,t}\\ {\mathbf{0}}_{1\times 3}& 1& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 1& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 1& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{ccccc}{\left({\mathbf{w}}_{0,t}^g\right)}_{\times }& {\mathbf{w}}_{0,t}^a& {\mathbf{0}}_{3\times 1}& {\mathbf{h}}_R\left({\tilde{\alpha }}_t\right){\mathbf{w}}_t^v& {\mathbf{h}}_R\left(\tilde{\alpha }\right){\mathbf{w}}_{i,t}^a\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 0\end{array}\right]\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \triangleq & f_{u_t}\left({\overline{X}}_t\right)-{\overline{X}}_t{\mathcal{L}}_{\mathfrak{g}}\left({\overline{\mathbf{w}}}_t\right)\hfill \end{array}$$

(13)

with ${\overline{\mathbf{w}}}_t:=\mathrm{vec}({\mathbf{w}}_{0,t}^g,{\mathbf{w}}_{0,t}^a,{\mathbf{0}}_{3\times 1},{\mathbf{h}}_R\left({\tilde{\alpha }}_t\right){\mathbf{w}}_t^v,{\mathbf{h}}_R\left(\tilde{\alpha }\right){\mathbf{w}}_{i,t}^a)$. The right-invariant error between two trajectories ${\overline{X}}_t$ and ${\widehat{\overline{X}}}_t$ is defined as

$${\eta }_t^r={\widehat{\overline{X}}}_t{\overline{X}}_t\triangleq \mathbf{exp}\left({\overline{\xi }}_t\right)\approx {\mathbf{I}}_d+{\mathcal{L}}_{\mathfrak{g}}\left({\overline{\xi }}_t\right).$$

(14)

where ${\mathcal{L}}_{\mathfrak{g}}:{\mathbb{R}}^{3N+9}\to \mathfrak{g}$ maps a vector to the corresponding element of the Lie algebra. The $f_{u_t}(·)$ can be shown to satisfy the group affine property. Therefore, the right-invariant error dynamics are trajectory independent and satisfy ${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\eta }_t^r=f_{u_t}\left({\eta }_t^r\right)-{\eta }_t^r{f}_{u_t}\left(I_d\right)+\left({\widehat{X}}_t{\mathcal{L}}_{\mathfrak{g}}\left({\mathbf{w}}_t\right){\widehat{X}}_t^{-1}\right){\eta }_t^r\triangleq g_{u_t}\left({\eta }_t^r\right)+{\mathcal{L}}_{\mathfrak{g}}\left({\widehat{\mathbf{w}}}_t\right){\eta }_t^r$ [18]. Considering the IMU biases,

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\widehat{\mathbf{d}}}_{i,t}-{\widehat{\mathbf{R}}}_t{\mathbf{R}}_t^{\top }{\mathbf{d}}_{i,t}\right)& \approx {\xi }_t^{v_i}+{\left({\widehat{\mathbf{d}}}_{i,t}\right)}_{\times }{\widehat{\mathbf{R}}}_t({\mathbf{w}}_{0,t}^g-{\zeta }_{0,t}^g),\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\widehat{\mathbf{v}}}_{i,t}-{\widehat{\mathbf{R}}}_t{\mathbf{R}}_t^{\top }{\mathbf{v}}_{i,t}\right)& \approx {\left(\mathbf{g}\right)}_{\times }{\xi }_t^R+{\left({\widehat{\mathbf{v}}}_{i,t}\right)}_{\times }{\widehat{\mathbf{R}}}_t({\mathbf{w}}_{0,t}^g-{\zeta }_{0,t}^g)\hfill \\ & +{\widehat{\mathbf{R}}}_t{\mathbf{h}}_R\left(\tilde{\alpha }\right)({\mathbf{w}}_{i,t}^a-{\zeta }_{i,t}^a).\hfill \end{array}$$

(16)

where ${\zeta }_{0,t}^g={\widehat{\mathbf{b}}}_{0,t}^g-{\mathbf{b}}_{0,t}^g$ and ${\zeta }_{i,t}^a={\widehat{\mathbf{b}}}_{i,t}^a-{\mathbf{b}}_{i,t}^a$. To let the invariant error satisfy the log-linear property [18], ${\overline{A}}_t$ is defined by $g_{u_t}\left(\exp \left(\overline{\xi }\right)\right)={\mathcal{L}}_{\mathfrak{g}}\left({\overline{A}}_t\xi \right)+\mathcal{O}(\parallel \overline{\xi }{\parallel }^2)$. Using the first-order approximation $\exp \left(\overline{\xi }\right)\approx I_d+{\mathcal{L}}_{\mathfrak{g}}\left({\overline{\xi }}_t\right)$, we can have

$$\overline{A}=\left[\begin{array}{ccccccccc}{\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -\widehat{\mathbf{R}}& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\\ {\left(\mathbf{g}\right)}_{\times }& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\left({\widehat{\mathbf{v}}}_0\right)}_{\times }\widehat{\mathbf{R}}& -\widehat{\mathbf{R}}& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\\ {\mathbf{0}}_3& I_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\left(\widehat{\mathbf{p}}\right)}_{\times }\widehat{\mathbf{R}}& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{I}}_3& -{\left(\widehat{\mathbf{d}}\right)}_{\times }\widehat{\mathbf{R}}& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\\ {\left(\mathbf{g}\right)}_{\times }& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& -{\left({\widehat{\mathbf{v}}}_i\right)}_{\times }\widehat{\mathbf{R}}& {\mathbf{0}}_3& -\widehat{\mathbf{R}}{\mathbf{h}}_R\left(\tilde{\alpha }\right)& {\mathbf{0}}_{3\times 9}\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 9}\end{array}\right]$$

(17)

with $\overline{\xi }={[{\xi }_R^{\top },{\xi }_v^{\top },{\xi }_p^{\top },{\xi }_{d_i}^{\top },{\xi }_{v_i}^{\top },{\zeta }_{b_0^g}^{\top },{\zeta }_{a_0^g}^{\top },{\zeta }_{a_i^g}^{\top }]}^{\top }$. The covariance matrix of the process noise for the augmented error dynamics is $\widehat{\overline{Q}}=\left[\begin{array}{cc}{\mathrm{Ad}}_{{\widehat{\overline{X}}}_t}& {\mathbf{0}}_{12\times 6}\\ {\mathbf{0}}_{6\times 12}& I_6\end{array}\right]\mathrm{Cov}\left({\overline{\mathbf{w}}}_t\right){\left[\begin{array}{cc}{\mathrm{Ad}}_{{\widehat{\overline{X}}}_t}& {\mathbf{0}}_{12\times 6}\\ {\mathbf{0}}_{6\times 12}& I_6\end{array}\right]}^{\top }$.

Using the measurements ${\tilde{\omega }}_i$ of the gyroscope on the foot i,

$${}^B{\mathbf{v}}_i={\mathbf{R}}_t^{\top }\left[{}^W{\omega }_i(\widehat{x},\tilde{\alpha }{,}^{C_i}{\omega }_i)\times \mathbf{s}(\widehat{x},\tilde{\alpha })\right],$$

(18)

where ${}^W{\omega }_i=\mathbf{R}\left({\widehat{\theta }}_k\right){\mathbf{h}}_R{\left(\tilde{\alpha }\right)}^{C_i}{\omega }_i$. ${}^B{\mathbf{v}}_i$ can be expressed in the right-invariant observation structure ${\mathbf{Y}}_t={\mathbf{X}}_t^{-1}\mathbf{b}+{\mathbf{V}}_t$ in [18] as follows

$$\left[\begin{array}{c}{}^B{\mathbf{v}}_i\\ 0\\ 0\\ 0\\ -1\end{array}\right]=\left[\begin{array}{ccccc}{\mathbf{R}}_t^{\top }& -{\mathbf{R}}_t^{\top }{\mathbf{v}}_0& -{\mathbf{R}}_t^{\top }\mathbf{p}& -{\mathbf{R}}_t^{\top }{\mathbf{d}}_i& -{\mathbf{R}}_t^{\top }{\mathbf{v}}_i\\ {\mathbf{0}}_{1\times 3}& 1& 0& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 1& 0& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 1& 0\\ {\mathbf{0}}_{1\times 3}& 0& 0& 0& 1\end{array}\right]\left[\begin{array}{c}{\mathbf{0}}_{3\times 1}\\ 0\\ 0\\ 0\\ -1\end{array}\right]+\left[\begin{array}{c}\widehat{\mathbf{R}}{\mathbf{w}}_{C_i}^v\\ 0\\ 0\\ 0\\ 0\end{array}\right].$$

Then,

$$\begin{array}{cc}\hfill \overline{Y}& ={\left[\begin{array}{ccccc}{h}_p^{\top }\left(\tilde{\alpha }\right)& 0& 1& -1& 0\\ \mathbf{J}\left(\alpha \right)\dot{\alpha }+{{(}^B\omega )}_{\times }{h}_p\left(\alpha \right)& 1& 0& 0& -1\\ {}^B{\mathbf{v}}_i^{\top }& 0& 0& 0& -1\end{array}\right]}^{\top },\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill \overline{H}& =\left[\begin{array}{cccccc}{\mathbf{0}}_3& {\mathbf{0}}_3& -I& I& {\mathbf{0}}_3& {\mathbf{0}}_{3\times 18}\\ {\mathbf{0}}_3& -I& {\mathbf{0}}_3& {\mathbf{0}}_3& I& {\mathbf{0}}_{3\times 18}\\ {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& {\mathbf{0}}_3& I& {\mathbf{0}}_{3\times 18}\end{array}\right],\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \overline{N}& =\mathrm{diag}\left\{\begin{array}{c}\widehat{\mathbf{R}}{\mathbf{J}}_{v_i}\mathrm{Cov}\left({\mathbf{w}}_{i,t}^v\right){\mathbf{J}}_{v_i}^{\top }\left(\tilde{\alpha }\right){\widehat{\mathbf{R}}}^{\top },\\ \widehat{\mathbf{R}}\mathrm{Cov}\left({\mathbf{w}}_i^f\right){\widehat{\mathbf{R}}}^{\top },\\ \widehat{\mathbf{R}}\mathrm{Cov}\left({\mathbf{w}}_{C_i}^v\right){\widehat{\mathbf{R}}}^{\top }\end{array}\right\}.\hfill \end{array}$$

(21)

where $\mathrm{diag}\{·,·\}$ denotes the block diagonal operation, ${\mathbf{w}}_i^f$ incorporates the uncertainty in encoder measurements, kinematic model, the effect of slip and gyroscope bias [19], and ${\mathbf{w}}_{C_i}^v$ represents the noise of the measured velocity for the i-th contact point.

4. Learning-Assisted Multi-IMU Proprioceptive State Estimation

A notable limitation of MIPO [10] arises from its reliance on yaw angle measurements provided by a motion capture system. This dependency arises because using motion capture data significantly enhances the accuracy of the state estimation model. However, solely relying on IMU measurements or EKF predictions results in poor accuracy due to the high noise levels in IMU data and the lack of a cross-referencing method to correct these measurements. This reliance on motion capture systems, which are typically unavailable in outdoor environments or general experimental setups, limits the model’s applicability. To address this limitation, our method employs a learning-based approach to predict yaw angles from available data sources, including IMUs, joint motor encoders, and foot contact sensors. This predictive capability enables the replacement of direct motion capture measurements, thus fitting seamlessly within the multi-IMU odometry framework and broadening the model’s applicability to outdoor and less controlled environments. Figure 1 shows the schematic of our proposed LMIPSE.

To train a model for predicting yaw angles online, we first collect a dataset denoted by $\mathcal{D}=\{{\mathbf{x}}^{\left(k\right)}=\{{\tilde{\omega }}_{i,t_j},{\tilde{\mathbf{a}}}_{i,t_j},{\tilde{\alpha }}_{i,t_j},{\tilde{\mathbf{f}}}_{i,t_j}|i=0,1,\cdots ,4;j=k,k-1,\cdots ,k-n_{ws}+1\},{\mathbf{y}}^{\left(k\right)}={\theta }_y^{\left(k\right)}|k=1,\cdots ,N_{ds}\}$ where ${\mathbf{x}}^{\left(k\right)}$ are the inputs of the model while ${\mathbf{y}}^{\left(k\right)}$ is the output the model learns to predict, ${\tilde{\mathbf{f}}}_{i,t_j}$ is the foot force sensor measurements, $n_{ws}$ is the time window size, and $N_{ds}$ is the data set size. The dataset $\mathcal{D}$ can be randomly split into a training set ${\mathcal{D}}_{train}$ and a testing set ${\mathcal{D}}_{test}$. It is noted that mocap is only needed for collecting the data of ${\theta }_y^{\left(k\right)}$ but is needless for state estimation. Many model architectures have been proposed for signal processing, and we consider a 1D version of the ResNet (ResNet1D) which has proven efficient in learning inertial odometry [14,20]. Then, the ResNet1D denoted by $f(\mathbf{x};w)$ with parameters w can be trained by minimizing the following Smooth L1 loss (aka Huber Loss)

$$L_{{\delta }_{th}}(\mathbf{y}-f(\mathbf{x};w))=\left\{\begin{array}{cc}{\displaystyle \frac{1}{2}}{(\mathbf{y}-f(\mathbf{x};w))}^2\hfill & \mathrm{for}|\mathbf{y}-f(\mathbf{x};w)|\le {\delta }_{th},\hfill \\ {\delta }_{th}\left(|\mathbf{y}-f(\mathbf{x};w)|-{\displaystyle \frac{{\delta }_{th}}{2}}\right)\hfill & \mathrm{for}|\mathbf{y}-f(\mathbf{x};w)|>{\delta }_{th}.\hfill \end{array}\right.$$

(22)

The Huber loss (22) is chosen for its robustness to outliers in data. It smoothly transitions to a linear loss for errors larger than a certain threshold ${\delta }_{th}$. This combination of the squared loss and the absolute loss allows the Huber loss to provide the benefits of both types of losses.

It is noted that NN design significantly affects the generalization of the data-driven models [21]. Besides directly predicting yaw angles, we consider using a NN denoted by $f^{\prime }(\mathbf{x};w^{\prime })$ to model the mismatch between the yaw data (${\theta }_y$) from the mocap system and the yaw estimates (${\widehat{\theta }}_{DK,y}$) by the body IMU dead reckoning (DK). To learn an accurate $f^{\prime }(\mathbf{x};w^{\prime })$, first, we transform the dataset $\mathcal{D}$ into ${\mathcal{D}}^{\prime }=\{{\mathbf{x}}^{\left(k\right)}=\{{\tilde{\omega }}_{i,t_j},{\tilde{\mathbf{a}}}_{i,t_j},{\tilde{\alpha }}_{i,t_j},{\tilde{\mathbf{f}}}_{i,t_j}|i=0,1,\cdots ,4;j=k,k-1,\cdots ,k-n_{ws}+1\},{\mathbf{y}}^{\prime \left(k\right)}=\delta {\theta }_y^{\left(k\right)}|k=1,\cdots ,N_{ds}\}$ where $\delta {\theta }_y^{\left(k\right)}={\theta }_y^{\left(k\right)}-{\widehat{\theta }}_{DK,y}^{\left(k\right)}$ with ${\widehat{\theta }}_{DK,y}$ denoting the estimated yaw angle by body IMU DK. Then, we optimize $w^{\prime }$ in $f^{\prime }$ by minimizing $L_{{\delta }_{th}}({\mathbf{y}}^{\prime }-f^{\prime }(\mathbf{x};w^{\prime }))$ on the transformed dataset ${\mathcal{D}}^{\prime }$.

After finding a model $f^{\prime ★}:=f^{\prime }({\mathbf{x}}^{\left(k\right)};w^{\prime ★})$ that achieves small testing errors on ${\mathcal{D}}_{test}^{\prime }$, we use the model $f^{\prime ★}$ to predict the correction factor $\delta {\theta }_y$ for the yaw angle online. The predicted correction factor is then added to the yaw estimate ${\widehat{\theta }}_{DK,y}$ of IMU DK to obtain the final corrected yaw angle prediction ${\widehat{\theta }}_y$. The corrected ${\widehat{\theta }}_y$ is used for the prediction step and update step of the EKF. Specifically, at time instant k, the prediction step is

$$\begin{array}{cc}\hfill & {\widehat{x}}_{k+1}^{-}=f({\widehat{x}}_k^{+},u_k),\hfill \end{array}$$

(23a)

$$\begin{array}{cc}\hfill & P_{k+1}^{-}=F_k{P}_k^{+}{F}_k^{\top }+G_k{Q}_k^E{G}_k^{\top },\hfill \end{array}$$

(23b)

where $F={\displaystyle \frac{\partial f^E}{\partial x}}$ with $f^E$ denoting the state transition function in (7), $G={\displaystyle \frac{\partial f}{\partial w^E}}$ with $w^E$ denoting the process noise in (7), and $Q^E$ is the covariance matrix of $w^E$; the update step is

$$\begin{array}{cc}\hfill & P_k^{+}={({\left(P_k^{-}\right)}^{-1}+H_k^{\top }{R}_k^{-1}{H}_k)}^{-1},\hfill \end{array}$$

(24a)

$$\begin{array}{cc}\hfill & {\widehat{x}}_k^{+}={\widehat{x}}_k^{-}+P_k^{+}{H}_k^{\top }{R}_k^{-1}\left({\overline{y}}_k-H_k{\widehat{x}}_k^{-}\right),\hfill \end{array}$$

(24b)

where $H={\displaystyle \frac{\partial h}{\partial x}}$, and ${\overline{y}}_k=y_k-h\left({\widehat{x}}_k^{-}\right)+H_k{\widehat{x}}_k^{-}$.

5. Experiments and Validation

This section validates the proposed learning-assisted multi-imu proprioceptive state estimation approach using simulation and real data of a legged robot.

5.1. Dataset Description

The dataset utilized for this study was collected using a Unitree Go 1 robot in a lab space equipped with a highly accurate mocap system, as detailed in the work [10]. The robot was equipped with several proprioceptive sensors, including a body MEMS IMU, twelve joint motor encoders, and four foot pressure contact sensors. An MPU9250 IMU was also mounted on each robot foot to capture additional inertial information.

During indoor trials, the robot operated on flat ground and moved at speeds ranging from 0.4 to 1.0 m/s. Sensor data, including linear accelerations and angular velocities from the IMUs, torque and angle measurements from the joint motor encoders, and orientation data and positional data from the motion capture system. To ensure consistency, all collected data were first resampled to the same frequency, aligned to the same start and end times, and normalized such that the starting location of the robot was at the origin. These data were stored in a rosbag file and are publicly available [10].

A key limitation of the dataset is its short duration of approximately 45 seconds and the simple flat terrain of the testing environment, which might not significantly affect the IMU data’s susceptibility to drift. To address these limitations and simulate a more practical scenario, Gaussian noise and a time-integrated Gaussian noise drift were added to the body IMU data. This setup was intended to assess the robustness of the state estimation algorithms under conditions that more closely mimic operational environments with potential IMU data degradation.

5.2. Experimental Setup

In our study, the mocap system’s pose measurements serve as the ground truth, providing a high-fidelity benchmark for evaluating the performance of various state estimation methods. The Multi-IMU Proprioceptive Odometry (MIPO), which integrates multiple IMUs mounted on the robot’s feet, serves as the baseline against which other methods are compared. Additionally, we explore the body IMU dead reckoning as a first baseline method, where angular velocity measurements are integrated over time to estimate orientation. However, due to IMU drift, this method typically yields less reliable results.

Impacts of NN Design

To address the limitations caused by dataset variability and enhance the model’s generalizability across different robotic movements, we implemented two CNN-based approaches based on a 1D ResNet architecture [22]. The first approach is the Multi-IMU CNN Angle Estimator (MI-CAE), which predicts yaw angles directly, using mocap-derived orientation measurements as ground truth. The second approach is the Multi-IMU CNN angle Correction Enhancer (MI-CCE), which aims to predict a correction factor for the IMU-derived yaw angles, with the ground truth being the discrepancy between the IMU dead reckoning and mocap measurements. The corrected yaw angles are then integrated into the MIPO framework to enhance state estimation. The efficacy of these approaches is subsequently evaluated by comparing the resulting position estimates against the mocap system data, thus ensuring a robust assessment of each approach’s ability to improve upon traditional IMU dead reckoning techniques in dynamic environments.

NN architecture: We used a ResNet1D as the backbone of the model for processing sensor data. The ResNet1D structure comprises 9 residual blocks, each containing a sequence of convolutional layers, batch normalization, and ReLU activation functions. The output of the model is configured to produce a single continuous value, which presents the predicted yaw angles or the correction factor for IMU yaw angle measurements. Table 1 summarizes the architecture of the ResNet1D model.

Training settings: From a training perspective, the model employs the Adam optimizer with an initial learning rate of 0.001 and a weight decay of 0.001 to prevent overfitting. Additionally, a ReduceLROnPlateau learning rate scheduler is integrated to adjust the learning rate based on the performance, specifically reducing the rate if no improvement in loss is observed over 10 epochs. This is complemented by using Huber Loss as the loss function. The entire model is trained for 50 epochs.

5.3. Results and Discussions

5.3.1. Yaw Prediction Model

Table 2 summarizes the testing accuracy of the two CNN-based models using different training set sizes. After experimenting with various sizes, we chose to use 3000 random data points to balance performance and generalization.

Furthermore, Figure 2 compares the yaw prediction results of MI-CAE and MI-CCE, demonstrating that MI-CCE achieved smaller prediction errors. Since the yaw predictions are used for EKF, we further tested the performance of the two models for state estimation.

5.3.2. Proprioceptive State Estimation

We assessed the performance of different state estimation methods under both low-noise and high-noise conditions. Low Noise Condition refers to the original IMU data measured using a high-quality IMU in a controlled lab environment. High Noise Condition is created by adding Gaussian noise and drift to the IMU measurements to mimic the unstable environment.

Table 3 summarizes the results of various state estimation methods under low IMU noise conditions. The MIPO method [10] achieves significantly higher positioning accuracy compared to SIPO. However, using yaw angles estimated by EKF rather than motion capture (mocap) measurements substantially increases the position estimation errors for both SIPO and MIPO. When using the body IMU dead reckoning (DK) method, the performance of MIPO with mocap is slightly degraded, and this degradation compounds over time. The MIPSE with RI-EKF reduces estimation errors but still results in much larger errors compared to MIPO using mocap. In contrast, the MI-CAE method uses CNN-predicted yaw angles derived from multiple IMU measurements for state estimation. This results in an average drift of 15.65%, a median drift of 15.87%, an RMSE of 0.304946, and a maximum RSE of 0.855203. The MI-CCE method improves upon MI-CAE by using a CNN to predict an angle correction factor from multiple IMU measurements, which is then used to correct the IMU data before performing state estimation. This method shows slightly lower drift percentages and error metrics, with an average drift of 15.65%, a median drift of 15.87%, an RMSE of 0.304946, and a maximum RSE of 0.855203. The discrepancy in state estimation performance is smaller than the difference in yaw prediction accuracy between MI-CAE and MI-CCE, demonstrating that prediction errors can be mitigated by the EKF.

Notably, while the MIPO method uses angle information from a motion capture system and achieves the lowest error metrics, the CNN-enhanced methods, MI-CAE and MI-CCE, exhibit similar levels of accuracy without the reliance on external motion capture systems.

Table 4 summarizes the results of various state estimation methods under high IMU noise. Our LMIPSE using CCE achieved the best results regarding all three performance measurements. Furthermore, LMIPSE is more robust to increased noise level than MIPO with DK, as the median drift was increased only by 0.03% for LMIPSE as compared to 0.25% for MIPO with DK. Additionally, Figure 3 compares the trajectory estimation, and Figure 4 compares the state estimation results of LMIPSE and other methods.

6. Concluding Remarks

A learning-assisted EKF approach was presented for multi-IMU proprioceptive state estimation of quadruped robots. In particular, a 1D CNN-based model for estimating the measurements required for EKF was trained offline using the ground-truth measurements from the mocap system and the buffers of proprioception data (including IMUs and joint encoders measurements). Then, the trained model was used to estimate the measurements online for EKF. Experiments demonstrated that the proposed approach improved the accuracy of the state estimation as compared to multi-IMU state estimation using RI-EKF and proprioceptive odometry without using multiple IMUs. In future work, we will consider learning-based approaches for tuning the hyperparameters of the filters [23].