Model-Based Slippage Estimation to Enhance Planetary Rover Localization with Wheel Odometry

Abstract

The exploration of planetary surfaces with unmanned wheeled vehicles will require sophisticated software for guidance, navigation and control. Future missions will be designed to study harsh environments that are characterized by rough terrains and extreme conditions. An accurate knowledge of the trajectory of planetary rovers is fundamental to accomplish the scientific goals of these missions. This paper presents a method to improve rover localization through the processing of wheel odometry (WO) and inertial measurement unit (IMU) data only. By accurately defining the dynamic model of both a rover’s wheels and the terrain, we provide a model-based estimate of the wheel slippage to correct the WO measurements. Numerical simulations are carried out to better understand the evolution of the rover’s trajectory across different terrain types and to determine the benefits of the proposed WO correction method.

1. Introduction

Robotic exploration of the surface of celestial bodies in the Solar System is a key phase of planetary science investigations. In situ measurements are fundamental to significantly expand our knowledge of surface formation and evolution, habitability, geology and geophysics. Rovers on the Martian surface have led to unprecedented measurements (e.g., [1]), studying the planet’s habitability and seeking signs of past microbial life [2,3].

To conduct science investigations across a wide region of the body’s surface, the navigation of a rover must fulfil tight engineering constraints. Safe mobility is one of the main requirements, and is directly related to the morphology of the explored area and the properties of the terrain [4]. Serious risks are caused by high-slippage terrains, slopes and hazards, which strongly affect locomotion. The NASA rover Spirit, for example, got stuck in soft soil leading to the descoping of the mission as a stationary science platform [5]. An accurate characterization of the soil properties is then fundamental to prevent these undesired motion failures. Onboard systems can then be designed to identify a terrain type by using images of the environment [6]. The predicted terrain class can then be used to predict the slippage conditions and correct nonsystematic errors of wheel odometry [7].

The reconstruction of the rover’s pose through the processing of encoder measurements is also limited by systematic errors [8,9]. WO and inertial measurement unit (IMU) data are directly used to compute a vehicle’s motion from wheel displacements. This approach is generally based on a simplified modeling of the kinematic equations, leading to a significant error propagation [10]. To enhance the determination of the rover’s pose additional measurements are collected by onboard instrumentation, as, for example, navigation cameras for visual odometry [11,12,13] and light detection and ranging (LiDAR) [14]. Images can also be used to estimate the rover’s slip by jointly measuring the pose of the wheels’ traces and the rate of turn [15].

A correction of WO and IMU measurements is fundamental to absorb mismodeling errors when imaging and ranging data are not available. A novel method was introduced by Ojeda et al. [16] to detect the wheel-slippage for WO corrections by measuring motor current. Terramechanics and regression functions were adopted by Yamauchi et al. [17] to improve the localization accuracy that is obtained using IMU-based systems. Censi et al. [18] developed a technique based on WO and an exteroceptive sensor to simultaneously calibrate odometry and sensor parameters. Here, we present a method to compute the velocity slippage correction, which is the difference between the linear speed of the vehicle and the relative speed of the wheel measured by WO [16], with an accurate modeling of a rover’s dynamics. To accurately determine the rover’s motion on different terrain types, we implemented a detailed wheel–soil interaction that accounts for rover’s wheel and terrain deformability [19,20,21]. Previous works were in general based on the assumption of a rigid wheel on different terrain types [10,22,23]. This hypothesis, however, is not well-suited for rovers equipped with flexible wheels, i.e., ESA ExoMars rover [24].

In this work, the rover was modeled by using size and mass of ESA ExoMars rover [25]. Different terrain types were included in our software, including cloddy, mixed drift-cloddy and drift soils. Their physical properties rely on Martian soil simulants retrieved from semi-empirical Earth-based observations [26,27]. Numerical simulations of different scenarios are reported to test the advantages of the proposed method for WO correction.

This paper is organized as follows: Section 2 provides an insight into the modeling of the rover navigation and defines the full set of equations of motion used to describe its dynamics; Section 3 presents the outcomes of the numerical simulations for rover navigation across different terrain conditions; Section 4 provides the description of the path planning strategy adopted to simulate navigation and assesses the accuracy of the WO correction method. Finally, Section 5 outlines the conclusions of the paper.

2. Methods

2.1. Dynamical Model and Slippage Estimation

The motion of a wheeled vehicle on unprepared ground, such as planetary surfaces, requires an accurate modeling of the compliant wheel on compliant ground. The resistance to motion encountered by the wheels is due to the energy dissipation in the wheels and the ground. Figure 1a shows a scheme of a compliant wheel that moves on a compliant ground. According to the model introduced in [19,20,28], we assumed that the deformation in the tire, $f$, occurs only if the maximum pressure exerted on the ground, ${\sigma }_z,$ exceeds the pressure, $p_{gr}$, needed to deform the tire. As represented in Figure 1a, this value is exceeded at point A of the contact zone.

We assumed that, in the first part of the contact, i.e., in the interval (${\theta }_1, {\theta }_0$] reported in Figure 1a, the wheel behaves as a rigid body and the pressure is linked to the sinkage. For $\theta \in$ [$-{\theta }_1, {\theta }_1$], we assumed a constant pressure equal to $p_{gr}$, which stands for an equivalent stiffness of flexible wheels. The normal stress, ${\sigma }_z,$ is defined as follows:

$${\sigma }_z=\{\begin{array}{c}\left(\frac{k_c}{b}+k_{\varphi }\right)h{\left(\theta \right)}^n for \theta \in \left({\theta }_1, {\theta }_0\right]\\ p_{gr} for \theta \in \left[-{\theta }_1, {\theta }_1\right]\end{array}$$

(1)

where $k_c$ and $k_{\varphi }$ are the cohesive and the frictional moduli that define the stiffness of the soil, $h\left(\theta \right)$ is the sinkage, b is the width of the wheel and $n$ is the exponent of terrain deformation that is $\le \frac{1}{2}$ and $\approx 1$ for cohesive and frictional soils, respectively.

The wheel–soil interaction results in a vertical force, $F_z$, a longitudinal force, $F_x$, a lateral force, $F_y$, and an aligning moment, $M_z$, acting on each vehicle’s wheel (Figure 1b). The vertical force $F_z$ is:

$$F_z=R_{w}b{{\displaystyle \int }}_{{\theta }_1}^{{\theta }_0}\left[{\sigma }_z\left(\theta \right)\cos \left(\theta \right)+{\tau }_x\left(\theta \right)\sin \left(\theta \right)\right] d\theta +2b{x}_A{p}_{gr}=\frac{mg}{n_w}\cos \left(i\right)$$

(2)

where $R_w$ is the radius of the wheel; ${\tau }_x$ is the longitudinal shear stress and $x_A$ is the x-coordinate of point A; m is the mass of the vehicle, $n_w$ is the number of the wheels, g is the gravitational acceleration and i is the inclination of the plane. The first contribution to $F_z$ in Equation (2) represents the vertical force due to the normal load, ${\sigma }_z$, directed radially along arc AB, which depends on the sinkage, $h\left(\theta \right)$, and to the vertical projection of the longitudinal shear stress ${\tau }_x$, acting in the contact area delimited by (${\theta }_1, {\theta }_0$]. The second contribution to $F_z$ is due to the constant pressure, $p_{gr}$, which acts in the z direction of the wheel reference frame {W} along line A’A. Solving the nonlinear Equation (2), one obtains the values of ${\theta }_0$ and ${\theta }_1$.

The longitudinal force ($F_x$) is:

$$F_x=R_{w}b{{\displaystyle \int }}_{{\theta }_1}^{{\theta }_0}[-{\sigma }_z\left(\theta \right)\sin \left(\theta \right)+{\tau }_x\left(\theta \right)\cos \left(\theta \right)]d\theta +R_{w}b{{\displaystyle \int }}_{-{\theta }_1}^{{\theta }_1}{\tau }_x\left(\theta \right)d\theta$$

(3)

where the first integral term of Equation (3) is the contribution of the radial pressure ${\sigma }_z$ and the longitudinal shear stress ${\tau }_x$ acting along arc AB and projected along the x direction of the wheel’s reference frame {W}; the second integral term is due to the action of ${\tau }_x$ along line A’A. The expression of ${\tau }_x$ at a generic point of the contact zone is assumed to be:

$${\tau }_x\left(\theta \right)=sgn\left(\sigma \right)\left[c+{\sigma }_z\tan \left(\varphi \right)\right]\left[1-e^{-\frac{R_w\left|\sigma \right|}{K_x}\left({\theta }_0-\theta \right)}\right]$$

(4)

and $c$, $K_x$, $\varphi$ are the cohesion, the modulus of the shear deformation in the $x$ direction and the internal friction angle, respectively. The longitudinal slip, $\sigma$, is a measure of the difference between the commanded speed of the vehicle, $\mathsf{\Omega }{R}_w$, controlled through the wheels’ angular speed, and its actual speed, $V$:

$$\sigma =\{\begin{array}{c}\frac{\mathsf{\Omega }{R}_w-V}{\mathsf{\Omega }{R}_w} in driving\\ \frac{\mathsf{\Omega }{R}_w-V}{V} in braking\end{array}$$

(5)

Lateral force is due to lateral slip, which is caused by a non-zero sideslip angle, $\alpha$, between the xz plane and the direction of the velocity of the center of the wheel. The lateral force ($F_y$) is expressed as:

$$F_y=-R_{w}b{{\displaystyle \int }}_{-{\theta }_1}^{{\theta }_0}{\tau }_y\left(\theta \right)d\theta -R_{w}b{{\displaystyle \int }}_{-{\theta }_1}^{{\theta }_0}{F}_{yb}\left[R_w-h\left(\theta \right)\cos \left(\theta \right)\right]d\theta$$

(6)

the first integral term of (6) represents the effect of the lateral shear stress, ${\tau }_y$, that acts along the contact area delimited by [−${\theta }_1, {\theta }_0$], and the second integral term yields the contribution of the bulldozing force per unit width $F_{yb}$:

$$F_{yb}=D_1\left[c h\left(\theta \right)+D_2\frac{\rho g{h}^2\left(\theta \right)}{2}\right]$$

(7)

where $\rho$ is the density of the soil and $D_1$ and $D_2$ are coefficients that depend on the friction angle $\varphi$. The lateral shear stress is given by:

$${\tau }_y\left(\theta \right)=sgn\left(\alpha \right)\left[c+{\sigma }_z\tan \left(\varphi \right)\right]\left[1-e^{-\frac{R_w}{K_y}\left(1-\left|\sigma \right|\right)\left({\theta }_0-\theta \right)\left|\tan \left(\alpha \right)\right|}\right]$$

(8)

where $K_y$ is the modulus of the shear deformation in the $y$ direction.

The aligning moment, which tends to align the midplane of the wheel with the direction of the velocity, is modeled as:

$$M_z=-R_w^{2}b{{\displaystyle \int }}_{-{\theta }_1}^{{\theta }_0}{\tau }_y\left(\theta \right)\sin \left(\theta \right)d\theta$$

(9)

Equations (4), (6), (8) and (9) are retrieved, accounting for a certain symmetry with respect to both slip and sideslip angles.

We implemented the so-called model for high-speed cornering [28,29] to a 3-axle car-like vehicle, whose properties are reported in

Table 1. Properties of the 3-axle simulated rover based on the ESA ExoMars mission [25]. $J_z$ is the rover’s moment of inertia, $x_{ij}$ is the distance between the centre of mass of the rover and the i-th axle, $y_{ij}$ is the distance from the longitudinal axis to the ij-th wheel, t is the distance between the kingpin axes of the wheels and $v_{max}$ and $v_c$ are the vehicle’s maximum and nominal speed.

Parameter [Unit]	Value
m [kg]	310
$J_z$ [kg ${\mathrm{m}}^2]$	93.0
${\mathrm{x}}_1$ [m]	0.76
${\mathrm{x}}_2$ [m]	0.14
${\mathrm{x}}_3$ [m]	−0.59
$y_{i1}$ [m]	0.6
$y_{i2}$ [m]	−0.6
L [m]	1.35
$R_w$ [cm]	14.25
b [cm]	12.0
t [m]	1.2
$p_{gr}$ [kPa]	14.6
${\mathrm{v}}_{max}$ [km/h]	0.1
$v_c$ [km/h]	0.04

.

The full set of equations of motion is:

$$\begin{gathered} \dot{u}=vr-g\sin \left(i\right)\sin \left(\mathsf{\psi }\right)+\frac{1}{m}{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^2\left[F_{x_{ij}}\cos \left({\mathsf{\delta }}_{ij}\right)-F_{y_{ij}}\sin \left({\mathsf{\delta }}_{ij}\right)\right] \\ \dot{v}=-ur-g\sin \left(i\right)\cos \left(\mathsf{\psi }\right)+\frac{1}{m}{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^2\left[F_{x_{ij}}\sin \left({\mathsf{\delta }}_{ij}\right)+F_{y_{ij}}\cos \left({\mathsf{\delta }}_{ij}\right)\right] \\ \dot{r}=\frac{1}{J_z}{\displaystyle \sum }_{i=1}^n{\displaystyle \sum }_{j=1}^2\{F_{x_{ij}}\left[-y_{ij}\cos \left({\mathsf{\delta }}_{ij}\right)+x_i\sin \left({\mathsf{\delta }}_{ij}\right)\right]+F_{y_{ij}}\left[y_{ij}\sin \left({\mathsf{\delta }}_{ij}\right)+x_i\cos \left({\mathsf{\delta }}_{ij}\right)\right]+M_{z_{ij}}\} \\ \dot{\mathsf{\psi }}=r \end{gathered}$$

(10)

where ${\delta }_{ij}$ is the steering angle of the ij-th wheel; $u$ and $v$ are the components of the rover’s velocity in the body reference frame, and $\dot{\psi }$ is the yaw rate that is angular velocity about the $z$ axis. The first two equations in Equation (10) provide the vehicle accelerations, $\dot{u}$ and $\dot{v}$, in the direction of x and y body axes, respectively; the contribution of the yaw rate, of the gravitational acceleration and the action of the longitudinal and lateral forces are also included. The third equation of Equation (10) yields the rover’s angular acceleration $\dot{r}$ produced by the action of the interaction forces and the aligning moments. To compute the rover’s trajectory two more equations must be introduced:

$$\left[\begin{array}{c}\dot{X}\\ \dot{Y}\end{array}\right]=\left[\begin{array}{cc}\cos \left(\psi \right)& -\sin \left(\psi \right)\\ \sin \left(\psi \right)& \cos \left(\psi \right)\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]$$

(11)

The slippage correction to correct the WO measurements is based on the combination of the motion Equations (10) and (11) and the interaction forces Equations (2), (3), (6) and (9). Our method consists in an iterative process that starts from the acquisition of the wheel’s angular rate and the vehicle yaw rate from WO and IMU, respectively. By initially assuming a zero slippage and a first guess on the rover’s state ($X, Y, \psi$), the rover’s state is propagated through iterative calculations of the slippage, the interaction forces and all the derivatives values shown in Equations (10) and (11). At each step, the effect of the soil properties, and, thus, of the interaction forces, impacts on the components of the rover’s linear and angular accelerations. This allows the estimation of the actual slippage with Equation (5) and improves the estimates of the rover’s speed and yaw rate.

2.2. Soil and Rover Properties

Each terrain type is associated with a specific set of parameters that characterize the wheel–ground interaction (i.e.,$c, \varphi , n, k_c,k_{\varphi }, K_x, K_y,$ and $\rho ,$ mentioned in Section 2.1). The numerical simulations presented in this work are carried out through a modeling of the Martian surface to mimic several operative conditions. Our mission scenarios account for cloddy, mixed drift-cloddy and drift soils, which are common terrain types on Mars. The physical properties of these soil are reported in

Table 2. Soil properties.

Parameter [Unit]	c [Pa]	ϕ [deg]	n [-]	kc [N/${\mathit{m}}^{\mathit{n}\mathbf{+}\mathbf{1}}$]	kϕ [N/${\mathit{m}}^{\mathit{n}\mathbf{+}\mathbf{2}}$]	Kx = Ky [M]	ρ [kg/m3]
Cloddy soil	170	37°	1	1400	820,000	0.016	1550
Mixed drift-cloddy soil	220	33.1°	1	1400	820,000	0.016	1350
Drift soil	530	26.4°	1	1400	820,000	0.016	1150

[26,27].

An accurate definition of the dynamical equations requires a 3D modeling of the rover. We implemented rover’s main body (Figure 2), wheels and joints by using geometrical properties of ESA ExoMars rover [25], as shown in Table 1. The rover has six wheels, each with its own steering mechanism. The control system of the actuators is supposed to realize the Ackermann condition [28,30].

3. Traversability Analyses

Numerical simulations of different mission scenarios are presented in this section to show the effects of the soil properties on the traversability. The motion path is reported in a local reference frame (X,Y) centered at the rover’s initial position. The X axis points eastwards, and the Y axis northwards. First, we compute the trajectory discrepancies that are retrieved by integrating the kinematic and dynamical equations. Mission scenarios are then simulated to better understand the impact of steep terrains on the rover’s traversability.

3.1. Trajectory Integration with the Dynamical Model

A first step in our analysis is a strict comparison between the rover’s trajectories obtained through the integration of the kinematic [31] and dynamical equations. The kinematic approach does not depend on the properties of the soil, which strongly impact on the rover’s dynamical forces. In this set of numerical simulations, we assumed a smooth flat cloddy soil. The rover is moving with a steering angle $\delta =2^{°}$ for the two front wheels. The commanded speed $\overline{V}$ is imposed to match the ExoMars nominal speed $v_c$ = 40 m/h [25]. The rover’s initial heading angle is ${\psi }_0=\pi /2$. By applying the dynamical model, an initial speed $u_0=v_c$ is adopted, whereas the rover’s speed is equal to $v_c$ for the kinematic model. The rover’s equations of motion were integrated for ~5.5 h to complete a loop with a radius $R\approx 38.66$ m.

Figure 3 shows the differences $\mathsf{∆}x, \mathsf{∆}y$ and $\mathsf{∆}\psi$ between the trajectory integrated with the kinematic and dynamical equations. The trajectory resulting from the kinematic model is along a precise circular path, since the interaction between the wheeled vehicle and terrain is neglected. Those differences are then interpretated as an error in the rover’s position. The kinematic and dynamical integrated trajectories differ by $\mathsf{∆}x=$ 7.89 m, $\mathsf{∆}y=10.20$ m and $\mathsf{∆}\psi ={17.80}^{°}$ at the end of the simulation. These significant errors in both position and orientation indicate the crucial role of an accurate dynamical modeling for rover’s navigation.

3.2. Traversability over Steep Terrains

The contribution of the terrain properties on the dynamical forces depends on the average slope of the surface. Steeper is the terrain larger are the path discrepancies for different soil types. By assuming constant terrain slopes and vehicle’s steering angles, we investigate the evolution of both uphill and downhill paths with different terrains.

A proportional derivative (PD) control scheme was implemented to control the wheels’ angular and linear velocity preventing from locking condition ($\overline{V}\approx v_{min})$ and exceeding the maximum speed ($\overline{V}\approx v_{max})$. If we assumed a constant commanded speed along the entire trajectory, the constraint on the maximum speed in the descending path would be violated, and the rover’s motion would result in a slow uphill traverse. Therefore, at each integration step the angular acceleration of the wheels, $\dot{\mathsf{\Omega }}$, and the value of the commanded velocity, $\overline{V}$, are given by:

$$\dot{\mathsf{\Omega }}=-K_p\left(V-V_d\right)-K_{d1}\dot{u}-K_{d2}\dot{v}$$

(12)

$$\overline{V}={\mathsf{\Omega }\mathrm{R}}_{\mathrm{w}}$$

(13)

which enable the fulfilment of the previous motion constraints by assuming $\dot{\mathsf{\Omega }}=0$ if $\overline{V}\ge v_{max}$ and $\dot{\mathsf{\Omega }}\ge 0$ or $\overline{V}\le v_{min}$ and $\dot{\mathsf{\Omega }}\le 0$. In our numerical simulations, the PD control gains $K_p$, $K_{d1}$ and $K_{d2}$ are 200, 50 and 50, respectively.

A first test case is based on a 10°-slope and a constant steering angle $\delta =4^{°}$ for the front and rear wheels. A stable rover speed $V_d=40$ m/h is obtained through extra drive torque when driving uphill. This improves the action of the tractive forces and increases the actual speed of the rover, almost maintaining the desired value. On the contrary, braking is applied during rover’s downhill path. In this case the commanded speed is reduced to obtain a negative value of the slip that induces a negative traction coefficient and a stronger braking action of the longitudinal forces.

Figure 4, Figure 5 and Figure 6 show the results of the trajectory integrated for 3 h by using an initial heading angle ${\psi }_0=\pi /2$ and speed $u_0=v_c$. At $t=0$ the commanded speed is set equal to the desired speed ($\overline{V}=V_d=40$ m/h). The rover’s speed is controlled within the boundary constraints for the three terrain types. The commanded speed approaches the value $v_{max}$ in the uphill phase, and tends to lower values when the rover goes downhill. The actual speed of the rover follows the desired value $V_d$. As expected, in the uphill the longitudinal slip is positive, while it tends to negative values in the downhill phases. The longitudinal slip reaches the maximum value for drift soil, $\sigma =0.153$, indicating high-slippage conditions.

A second case is based on the trajectory integration with the same vehicle’s steering angle and initial state conditions on a steeper surface with a 20° slope. This terrain inclination is an upper limit constraint on ExoMars traverse. Its onboard path planning standard will avoid slopes $\ge {20}^{°}$ [25]. The rover’s paths with a steeper terrain are dramatically different compared to the previous scenario (Figure 7). The rover completes one and a half loop on cloddy soil but its trajectory is along a straight line for mixed drift cloddy and drift soils. A steering angle of 4°, therefore, is not sufficient to follow a circular path because of different terrain parameters that affect the wheel–soil interaction. Figure 8 shows that the actual speed of the rover does not violate neither the maximum velocity nor locking constraints for the three terrain types. The longitudinal slip is not larger than 0.56, even in the case of drift soil (Figure 9). Uphill and downhill paths are then feasible for the rover over the simulated terrain models with a slope angle lower than 20°.

By simulating further terrain slopes, we identified the trafficability limit for each terrain type described in Table 2. This parameter is defined as the maximum slope angle that leads to a slip ratio close to unity [32]. Our results suggest that the trafficability limits are: (1) $~{30}^{°}$ for cloddy soil; (2) slope $~{26}^{°}$ for mixed drift-cloddy soil; and $~{21}^{°}$ for drift soil. These estimates are fully consistent with ExoMars expected performances [25].

4. Wheel Odometry Correction for Rover’s Localization

4.1. Simulated Trajectory

The dynamical model was used to define the trajectory of the rover in a Martian-like environment characterized by drift soil and hazards randomly distributed on the ground. The rover task is to reach the target point $\left({\mathrm{X}}_{\mathrm{f}}^{*}, Y_f^{*}\right)=\left(100.0, 350.0\right)$ m in the local reference frame avoiding obstacles along the path. We assumed to perfectly know the location of each object on the ground. To generate the rover’s reference trajectory for the numerical simulations on the localization method with WO correction, we implemented a simplified PD control that is a simplified approach compared to other strategies (e.g., [33,34,35]).

The initial condition of the rover’s trajectory consists of heading angle ${\psi }_0=0^{°}$, commanded speed $\overline{V}=$ 40 m/h and actual speed V =$\overline{V}$. To reach the goal, we used the PD control on the commanded speed presented in Section 3.2, and a control with variable steering strategies [30]. The integration of the dynamical equations allows us to determine the longitudinal slip and the value of the angle ${\theta }_0$. The values of the interaction forces are evaluated at each step accordingly to Equations (3), (6) and (9).

The first trajectory leg $(t\le 4.2$ h) is a straight line that is followed by tuning the steering angle as follows:

$$\delta =-K_{d}d+K_h\left({\psi }^{*}-\psi \right)$$

(14)

where d is the distance between the rover position and the desired path, ${\psi }^{*}$ is the desired heading angle and $K_d$ and $K_h$ are the control gains equal to 0.08 and 0.02, respectively. The desired speed is $V_d=0.4 v_{max}$. Then, this speed is kept constant with a steering angle $\delta =1^{°}$ for ~3.3 h. The last segment of the trajectory only changed the steering angle to $0^{°}$. Once the distance between the rover’s position and the target is less than 20 m, a proportional control on the steering angle and on the desired speed is adopted to reduce $V_d$ in approaching the target, which is reached after 13.2 h.

Figure 10 shows the path followed by the rover to enable hazard avoidance. The resulting trajectory is significantly longer than the straight line connecting the starting and target positions (shortest path). The rover’s speed is shown in Figure 11, together with the commanded and maximum speed. The evolution of the steering angle ${\delta }_{1R}$ of wheel 1R is reported in Figure 12.

4.2. Rover’s Localization

The reference trajectory was then used to run numerical simulations for the localization system presented in this study. The sensors that are modeled in our software provide WO encoder’s measurements of the wheels’ angular rate and IMU measurements of the heading direction and the yaw rate [36,37]. This dead-reckoning system is affected by non-systematic effects on demanding terrains since WO yields the travelled distance without considering the slippage between the wheels and the ground. Systematic errors due to the integration of these data in a simplified kinematic model accumulate in time leading to poor reconstructed trajectories. By applying our dynamical model presented in Section 2, we compute a slippage correction for WO to enhance the rover’s localization with dead-reckoning measurements only.

To determine the errors induced by WO and IMU data, we implemented a dead-reckoning estimator based on a kinematic model. The noise models for WO and IMU measurements are zero mean Gaussian. WO data are perturbed with a variable standard deviation ${\sigma }_{\omega }$ that increases with the wheel’s angular rate. A maximum standard deviation value of $0.25$ cm/s was adopted for a wheel’s angular rate of 0.0779 rad/s (i.e., $\overline{V}=\omega R_w=1.11$ cm/s). The IMU yaw rate measurements are affected by ${\sigma }_{\dot{\psi }}=0.013$ deg/s. Both instruments are assumed to collect data with a 60-s sampling rate. Errors associated with the terrain properties and the simplified motion equations lead to significant discrepancies between the simulated and reconstructed trajectory with uncorrected WO+IMU observations. The final estimated position of the rover is ~60 m away from the target location, confirming the poor performance of this dead-reckoning system.

Our localization system is based on the modeling of the dynamical equations to significantly reduce systematic effects. The wheels angular rate $\tilde{\omega }$ from WO and the vehicle’s yaw rate $\tilde{\dot{\psi }}$ from IMU enable an update of the rover’s state $\tilde{q}={\left[\tilde{X}, \tilde{Y}, \tilde{\psi }, \tilde{\dot{\psi }}, \tilde{u}, \tilde{v}\right]}^T$ through Equations (10) and (11). The initial rover’s coordinates and heading angle are assumed to be perfectly known as in the previous case for uncorrected WO+IMU. The initial components of the rover’s speed are assumed to be ${\tilde{u}}_0={\tilde{\omega }}_0{R}_w$ and ${\tilde{v}}_0=0$. At each time step, $\mathsf{\Delta }t=60$ s, the measured value $\tilde{\omega }$ and $\tilde{\dot{\psi }}$ are integrated in the dynamical equations to enable a slippage correction. The estimated $\tilde{X}$, $\tilde{Y}$, $\tilde{\psi }$, $\tilde{u}$ and $\tilde{v}$ are then propagated for $60$ s until new measurements of WO and IMU are available.

To account for mismodeling in the terrain properties that significantly affect the wheel–soil interaction, we introduced a 10% error in the soil parameters. We assumed to detect drift soil with cohesion, $c=583$ Pa, friction angle, $\varphi ={23.76}^{°}$, and density, $\rho =1035$ kg/m³. A major error contribution is related to our knowledge of the slope angle over the site. We then perturbed the slope in our dynamical model by including a variable slope angle between 0.1° and −0.1° and a white Gaussian noise with a standard deviation of ${\sigma }_i={0.1}^{°}$.

The slippage correction method provides an accurate estimation of the rover’s trajectory (Figure 13). This method allows us to compensate WO and IMU errors. Furthermore, the use of the dynamical equations reduces systematic errors leading to better estimates of the vehicle speed and yaw rate. The error in the final rover’s position is three times lower than uncorrected WO measurements (Figure 14), although we introduced errors in the terrain parameters and slope (Figure 15). Figure 16 shows the simulated vehicle’s speed and yaw rate in comparison with the estimates produced by WO correction method.

5. Summary

In this paper, we presented a method to determine the wheels’ slippage through an accurate modeling of the rover’s dynamical equations. This technique enables a model-based correction of WO measurements that are significantly affected by non-systematic effect associated with the terrain properties. An accurate model of the soil is fundamental to enhance the navigation of rovers across off-road terrains. The interaction forces between the rover’s wheels and the soil are used to predict the slippage condition that affect the accuracies of WO data. We modeled these forces by assuming compliant wheel on compliant ground.

Numerical simulations of different mission scenarios were carried out by considering a 3D model of the rover that operates on a Martian-like surface. We analyzed the rover’s trajectory evolution by accounting for different terrain types and conditions. Homogeneous and constant soil characteristics are assumed across the rover’s path. Our results are consistent with ExoMars expectations in terms of traversability of steep terrains, showing 21° and 26° slope limits for very fine sand terrain and coarse sand, respectively.

The dynamical equations are integrated in a localization system based on the processing of WO and IMU data only. We presented a method to correct the WO measurements by predicting the wheels’ slippage. This technique enables significant improvements in the estimates of the WO-based localization system. By perturbing the terrain properties in the dynamical model, we demonstrated that our method is robust to a partial lack of knowledge regarding soil physical parameters and slope angle. This method only requires a rough a priori characterization of the environment that could also be based on measurements collected by orbiters.