Online terrain estimation for autonomous vehicles on deformable terrains

TL;DR

This paper presents a real-time terrain estimation framework for autonomous vehicles on deformable terrains, using a reduced-order nonlinear terramechanics model and an unscented Kalman filter to improve navigation accuracy.

Contribution

It introduces a computationally efficient surrogate terramechanics model extending the Bekker model and demonstrates its effectiveness in real-time sinkage parameter estimation.

Findings

Estimated sinkage exponent within 4% of true value for clay and sandy loam.

Reduced prediction errors of vehicle states by orders of magnitude using the estimated parameters.

Surrogate model reduces computation cost by an order of magnitude compared to traditional SCM.

Abstract

In this work, a terrain estimation framework is developed for autonomous vehicles operating on deformable terrains. Previous work in this area usually relies on steady state tire operation, linearized classical terramechanics models, or on computationally expensive algorithms that are not suitable for real-time estimation. To address these shortcomings, this work develops a reduced-order nonlinear terramechanics model as a surrogate of the Soil Contact Model (SCM) through extending a state-of-the-art Bekker model to account for additional dynamic effects. It is shown that this reduced-order surrogate model is able to accurately replicate the forces predicted by the SCM while reducing the computation cost by an order of magnitude. This surrogate model is then utilized in a unscented Kalman filter to estimate the sinkage exponent. Simulations suggest this parameter can be estimated within 4% of its true value for clay and sandy loam terrains. It is also shown that utilizing this estimated parameter can reduce the prediction errors of the future vehicle states by orders of magnitude, which could assist with achieving more robust model-predictive autonomous navigation strategies.

Figures

Tables

Table 1. TABLE I : SCM terrain parameters.

Parameter	Symbol	Unit
Cohesive modulus	$k_c$	$\text{N}/{\text{m}}^{n+1}$
Frictional modulus	$k_{\varphi }$	$\text{N}/{\text{m}}^{n+2}$
Sinkage exponent	$n$	-
Shear deformation modulus	$k$	m
Cohesion	$c$	Pa
Angle of internal friction	$\varphi$	rad

Table 2. TABLE II : Terrain parameters for sand [ 8 ] .

Parameter	Value
$k_c$	1000 ($\text{N}/{\text{m}}^{n+1}$)
$k_{\varphi }$	1528600 ($\text{N}/{\text{m}}^{n+2}$)
$n$	1.08 (–)
$k$	0.024 (m)
$c$	200 (Pa)
$\varphi$	0.4712 (rad)

Table 3. TABLE III : Wheel states for benchmark simulation.

State	Value
Normal load	2500 (N)
Longitudinal slip	0.2 (–)
Camber angle	0 (rad)
Speed	5.5 (m/s)

Table 4. TABLE V : Measurement standard deviations used for sensor simulation.

State	Noise ($\sigma$)
$x$	1.2 (m)
$y$	1.2 (m)
$\psi$	0.0175 (rad)
$u$	0.25 (m/s)
$v$	0.25 (m/s)
${\omega }_z$	0.0175 (rad/s)

Table 5. TABLE VI : Terrain parameters for simulated terrains [ 7 ] .

Parameter	Sandy Loam	Clay
$k_c$	5300 ($\text{N}/{\text{m}}^{n+1}$)	13200 ($\text{N}/{\text{m}}^{n+1}$)
$k_{\varphi }$	1515000 ($\text{N}/{\text{m}}^{n+2}$)	692200 ($\text{N}/{\text{m}}^{n+2}$)
$n$	0.7 (–)	0.5 (–)
$k$	0.025 (m)	0.01 (m)
$c$	1700 (Pa)	4140 (Pa)
$\varphi$	0.5061 (rad)	0.2269 (rad)

Table 6. TABLE VII : Initial guess, estimated value, and estimation errors of the sinkage exponent n 𝑛 n for simulated terrains.

Terrain	Initial guess	True val.	Converged val.	% error
Sandy loam	0.9	0.7	0.722	3.1%
Clay	0.7	0.5	0.519	3.8%

Table 7. TABLE VIII : Mean squared error over entire simulation with varying prediction horizons for clay using estimated terrain parameter ( n = 0.519 𝑛 0.519 n=0.519 ) and initial guess ( n = 0.7 𝑛 0.7 n=0.7 ).

Time horizon	0.5 (s)		2.5 (s)		5 (s)
State	n=0.519	n=0.7	n=0.519	n=0.7	n=0.519	n=0.7
$x$	0.0034 (m)	0.0025 (m)	0.037 (m)	0.01 (m)	0.075 (m)	0.078 (m)
$y$	0.0035 (m)	0.0051 (m)	0.022 (m)	0.15 (m)	0.12 (m)	0.34 (m)
$\psi$	3.17e-05 (rad)	1.8e-04 (rad)	2.45e-04 (rad)	0.0089 (rad)	9.9e-04 (rad)	0.0115 (rad)
$u$	6.46e-05(m/s)	6.46e-05 (m/s)	1.33e-04 (m/s)	1.33e-04 (m/s)	1.2e-04(m/s)	1.2e-04 (m/s)
$v$	0.0027 (m/s)	0.013 (m/s)	0.0047(m/s)	0.15 (m/s)	0.005(m/s)	0.28 (m/s)
${\omega }_z$	6.2e-04 (rad/s)	0.004 (rad/s)	9.01e-04 (rad/s)	0.023 (rad/s)	9.05e-04 (rad/s)	0.04 (rad/s)

Equations

$$\sigma =(k_c/b+k_{\varphi })h^n$$

$$\tau ={\tau }_{\text{max}}(1-e^{-j/k})$$

$${\tau }_{\text{max}}=(c+\sigma \tan \varphi )$$

$$\sigma (\theta )=(k_c/b+k_{\varphi })h(\theta {)}^n$$

$$\tau (\theta )={\tau }_{\text{max}}(1-e^{-j(\theta )/k})$$

$${\tau }_{\text{max}}=(c+\sigma (\theta )\tan \varphi )$$

$$h(\theta )=\{\begin{array}{ll}r(\cos \theta -\cos {\theta }_{\text{f}})& {\theta }_{\text{m}}\le \theta \le {\theta }_{\text{f}}\\ r(\cos {\theta }_{\text{e}}-\cos {\theta }_{\text{f}})& {\theta }_{\text{r}}\le \theta \le {\theta }_{\text{m}}\end{array}$$

$${\theta }_{\text{m}}=(a_0+a_{1}s){\theta }_{\text{f}}$$

$${\theta }_{\text{f}}={\cos }^{-1}(1-h_f/r)$$

$${\theta }_{\text{e}}={\theta }_{\text{f}}-(\theta -{\theta }_{\text{r}})({\theta }_{\text{f}}-{\theta }_{\text{m}})/({\theta }_{\text{m}}-{\theta }_{\text{r}})$$

$${\theta }_{\text{r}}={\cos }^{-1}(1-\Lambda h/r)$$

$$j(\theta )=\{\begin{array}{ll}r[({\theta }_{\text{f}}-\theta )-(1-s)(\sin {\theta }_{\text{f}}-\sin \theta )]& s\ge 0\\ r[({\theta }_{\text{f}}-\theta )-(1/(1+s))(\sin {\theta }_{\text{f}}-\sin \theta )]& s<0\end{array}$$

$$h_0={\left[\frac{3W}{b(3-n)(k_c/b+k_{\varphi })\sqrt{2r}}\right]}^{\frac{2}{2n+1}}$$

$$F_z={\int }_{{\theta }_{\text{r}}}^{{\theta }_{\text{f}}}rb(\tau (\theta )\sin (\theta )+\sigma (\theta )\cos (\theta ))d\theta$$

$$h_0^{\prime }=h_0-F_z(h_0)/F_z^{\prime }(h_0)$$

$$F_y={\int }_{{\theta }_{\text{r}}}^{{\theta }_{\text{f}}}rb{\tau }_y(\theta )$$

$${\tau }_y(\theta )={\tau }_{\text{max}}(1-e^{-|j_y(\theta )|/k_y})$$

$$j_y(\theta )=r(1-s)({\theta }_{\text{f}}-\theta )\tan \beta$$

$${\tau }_y(\theta )={\tau }_{\text{max}}(j/k_y)(e^{1-j_y(\theta )/k_y})$$

$${\tau }_y^{\ast }(\theta )={\tau }_{\text{max}}(1-e^{-|j_y^{\ast }(\theta )|/k_y^{\ast }})g_1(v,s,F_z,n)$$

$${\tau }_y^{\ast }(\theta )={\tau }_{\text{max}}(j/k_y^{\ast })(e^{1-j_y^{\ast }(\theta )/k_y^{\ast }})g_1(v,s,F_z,n)$$

$$\displaystyle\begin{split}j^{*}_{y}(\theta)=&-|r(1-s)(\theta_{f}-\theta)\tan\beta|+\\ &\text{sign}(\beta)\big{(}r\sin(\theta)\Delta\delta\,g_{2}(v,s,F_{z},n)\big{)}\end{split}$$

$$k_y^{\ast }=g_3(v,s,F_z,n)$$

$$\dot{z_b}=\left[\begin{array}{c}u\cos \psi -(v+L_{\text{f}}{\omega }_z)\sin \psi \\ u\sin \psi +(v+L_{\text{f}}{\omega }_z)\cos \psi \\ w_z\\ a_x\\ (F_{\text{yf}}+F_{\text{yr}})/M_t-u{\omega }_z\\ (F_{\text{yf}}{L}_{\text{f}}-F_{\text{yr}}{L}_{\text{r}})/I_{zz}\end{array}\right]$$

$$z_b:=\left[\begin{array}{c}x\\ y\\ \psi \\ u\\ v\\ {\omega }_z\end{array}\right]=\left[\begin{array}{c}\text{global x position of front axle}\\ \text{global y position of front axle}\\ \text{yaw angle}\\ \text{longitudinal velocity}\\ \text{lateral velocity }\\ \text{yaw rate}\end{array}\right]$$

$$z:=\left[\begin{array}{c}{z}_b\\ n\end{array}\right],\dot{z}=\left[\begin{array}{c}\dot{z_b}\\ 0\end{array}\right]$$

$$z_{k+1}=F(z_k,v_k)$$

$$y_k=H(z_k,n_k)$$

$$Z_{k-1}=[\hat{z}\hat{z}\pm (\sqrt{(L+\lambda )P_z}{)}_i]$$

$$\lambda ={\alpha }^2(L+\kappa )-L$$

Full text

I Introduction

Autonomous ground vehicles (AGVs) have drawn interest for military applications to perform tasks, such as supply transport, in unsafe environments that could pose a threat to human operators [1]. Three considerations about military AGVs are important to motivate this work. First, military vehicles often need to operate off-road on deformable terrains, where the vehicle’s mobility is dependent on the highly nonlinear tire forces generated at the tire-terrain interface [2]. Second, increasing the mobility of military AGVs is a critical need [3]. Third, state-of-the-art approaches to navigate such vehicles typically rely on model dependent architectures, such as Model Predictive Control (MPC) [3, 4]. Therefore, when the AGVs are operated on deformable terrains, a more accurate knowledge of the terrain parameters becomes a critical enabler to maximize the mobility of the AGVs.

Much research has been performed in developing terramechanics models for off-road applications, which can be divided into empirical models, physics-based models, semi-empirical models [2]. Empirical model are the simplest; however, such models do not generalize well beyond the experimental test conditions used for their development. On the other hand, physics-based finite and discrete element models have proven to be of the highest fidelity, but the large computational efforts required renders them infeasible for real-time tire force prediction, thus limiting their applicability for use in AGVs and real-time terrain estimation [2]. More promising candidates, and perhaps the most widely used, for real-time tire force prediction on deformable terrains are the semi-empirical models based upon the classical terramechanics theory developed by Bekker, including the Soil Contact Model [5, 6, 7, 8]. In these models, the tire is typically assumed rigid and the deformation is assumed to occur only in the terrain [7]. To model the complex tire-terrain interactions, these terramechanics models rely on knowledge of terrain-specific parameters such as cohesion, internal friction angle, or sinkage exponent. During vehicle operation, these parameters may not be explicitly known or may be varying due to non-uniform terrains. Therefore, real-time terrain estimation is necessary in AGVs to improve the accuracy of the terramechanics models online and generate better informed control commands. Having this capability would also provide insight into traversability of terrains, such that path planning algorithms can reroute the vehicle to avoid regions where loss of mobility or excessive power consumption is likely to occur [9].

Researchers have already recognized this need and a limited number of results are available in the literature. In particular, in [5, 10], a Bayesian procedure is utilized for terrain parameter identification, but making this approach work online is subject to future research. Other researchers have proposed an online algorithm for estimating soil cohesion and internal friction angle utilizing a linear least-squares estimator for a rover [1, 11]. The algorithm relies on simplifying classical terramechanics equations through linear approximations to increase computational efficiency and subjects the rover to periodic high and low speed traverses [1]. However, linear approximations can lead to inaccurate stress approximations [12], and hence inaccurate force prediction, and periodically operating at low speeds is not desirable when maximum mobility is desired. Hence, online estimation of deformable terrain parameters for off-road AGVs is still an open research area and is the focus of this work.

This study presents a new approach for online terrain parameter estimation. First, due to the large computation time associated with integrating stresses in SCM and limitations of classical terramechanics equations, a nonlinear reduced-order model is developed by extending the work presented in [12] to account for additional dynamic effects such that a sufficient agreement with SCM can be achieved. Then, the reduced-order terramechanics model is incorporated in a 3 DoF bicycle model [4] to create an estimation model, whereas the actual vehicle is represented with a 11 DoF plant model with SCM. The predictions from the estimation model are fused with measurements from the plant model in an Unscented Kalman Filter (UKF) to identify the dominant terrain parameter, namely, the sinkage exponent. The result is an online terrain estimation approach that can be used to better inform control and path-planning algorithms for AGVs.

The rest of this paper is organized as follows. Sec. II first briefly reviews the SCM model used in the plant model. Then a state-of-the-art fast terramechanics model used as a benchmark is introduced and the significant deviations of its predictions from SCM are demonstrated. This model is then modified to improve its accuracy vis-à-vis SCM, so that a suitable estimation model is obtained. Sec. III presents the vehicle models, both the plant model as well as the estimation model. The terrain estimation procedure based on UKF is summarized in Sec. IV. Sec. V gives the simulation results including the accuracy of the estimations and their ability to improve the predictive accuracy of the 3 DoF model. Finally, Sec. VI gives the conclusions drawn from this work.

II Terramechanics Models

II-A Soil Contact Model (SCM)

This section briefly reviews the terramechanics model adopted in this work to represent the tire-terrain interactions in the plant simulations with high fidelity. This model is also used to evaluate the accuracy of the fast terramechanics models, including a state-of-the-art model and the surrogate model developed in this work. As such, this model serves as the ground truth for the purposes of this work.

The terramechanics model used in this study for generating the lateral tire forces acting on the vehicle is based on the Soil Contact Model (SCM) reported in [5, 13]. Verification of the model can be found in [14]. The SCM calculates relevant forces and torques acting on a 3 dimensional object in contact with a deformable terrain as summarized below.

The SCM algorithm relies on a discretized mesh of the tire and terrain to search for contact points at the tire-terrain interface. In the contact detection step, the vertices of the tire mesh are projected onto the nearest vertices of the terrain digital elevation map, effectively arranging the contact vertices in individual columns. The sinkage at each vertex can be determined from the minima of each column, assuming the vertex location is a point of sinkage. The effective contact width, $b$, can then be determined from the footprint’s area and contour length [13].

Following contact detection, the algorithm calculates the stresses at each contact node of the footprint as follows. The pressure, $\sigma$, is expressed as [15]

[TABLE]

The shear stress, $\tau$, is expressed as [16]

[TABLE]

with $\tau_{\text{max}}$ given as

[TABLE]

In the above expressions $h$ is the sinkage, $b$ is the tire effective width, and $j$ is the shear deformation. The remaining parameters are internal parameters characterizing the terrain as summarized in Table I. The forces generated at the tire-terrain interface can then be given by integrating the stresses over the entire contact patch. The above overview is a summary of [5]; a more complete discussion is given in [13].

SCM is a rather complex model due to the discretizations and integrations involved and may thus not be suitable for real-time parameter identification purposes. It has been shown that the accuracy of SCM is heavily influenced by the discretization resolution [13]. Furthermore, several SCM operations are of $N^{2}$ complexity, where $N$ is the number of grid nodes [13]. As an example, for a discretization of just 200 total nodes (100 per tire in a bicycle model), the time required by SCM can be around 20 ms [13]. Furthermore, for the UKF, the estimation method used in this work, 17 sigma points must be generated as discussed in Sec. IV, each calling the terramechanics model twice (once per tire in the bicycle model). Thus the total time spent calculating tire forces can be around 350 ms per a single UKF iteration. Finally, taking into account that many UKF iterations are needed to achieve estimation convergence, a UKF with SCM can be expected to take several minutes to converge, which is impractically long. Therefore, faster terramechanics models are needed.

II-B State-of-the-Art Fast Terramechanics Model

Much less computationally demanding solutions better suited for online estimation are given by Bekker-based models. These models are again based on (1)-(3); however, $\sigma$, $\tau$, $h$, and, $j$ are now replaced by functions of the angle of contact, $\theta$. As such, (1)-(3) are rewritten as:

[TABLE]

[TABLE]

[TABLE]

where $h(\theta)$ is given by

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

[TABLE]

where $r$ is the radius of the tire; $\theta_{\text{f}}$ is the angle at which the front of the tire comes into contact with the terrain; $\theta_{\text{m}}$ is the location of maximum normal stress with $a_{0}$ and $a_{1}$ as terrain parameters typically taking on values of 0.4 and 0–0.3, respectively [17]; $s$ is the longitudinal slip of the tire; $\theta_{e}$ is the equivalent front contact angle for angles less than $\theta_{m}$; $\theta_{\text{r}}$ is the angle at which the rear of the tire loses contact with the terrain; and $\Lambda$ is a property of the terrain characterizing the sinkage ratio.

Finally, $j(\theta)$ is given as

[TABLE]

The maximum sinkage can be calculated in an iterative fashion by using the Newton-Raphson method as proposed in [8] as follows. The maximum sinkage is initialized as the static sinkage, which is based on the load on the tire $W$:

[TABLE]

However, due to dynamic effects, such as slippage, additional sinkage is induced. To account for this, the reaction force is calculated as

[TABLE]

and a new sinkage is determined as

[TABLE]

The iterative procedure terminates when the calculated reaction force is within a specified tolerance of the normal force applied to the tire. Once the maximum sinkage is determined, the lateral force $F_{y}$ can then be calculated in a similar fashion as in [6, 8], i.e.,

[TABLE]

with $\tau_{y}(\theta)$ given as

[TABLE]

where

[TABLE]

and $\beta$ is the side slip angle.

Depending on the soil type $\tau_{y}(\theta)$ can also be represented with different a formulation such as

[TABLE]

Other formulations can be found in [7]. All relevant variables are depicted in Fig. 1.

To assess the accuracy of (16) compared to SCM, a simulation is run in Chrono [18]. The simulation utilizes a single wheel test bed operating on a sand-like terrain using Chrono’s built-in SCM terrain. The test bed allows for individual control of the tire’s velocity, load, longitudinal slip, and lateral slip. The terrain properties used in this simulation are representative of sand and given in Table II. The simulation sweeps the tire through a range of lateral slips with a 1 Hz sine wave. The load, longitudinal slip, camber angle, and linear velocity of the tire are all held at the constant values given in Table III.

Fig. 2 shows the results of an SCM simulation run in Chrono (orange) and the force predicted by (16) (blue). The term $k$ is assumed to be constant, rather than a function of lateral slip as in [6]. This is to maintain consistency with the SCM formulation used in Chrono. As seen in the figure, the base model of (16) captures the overall trend, at least in the linear region around zero lateral slip, but averages the two distinct curves seen in SCM. This is because the current formulation does not account for the hysteresis effects of varying lateral slip; i.e., the shear deformation of (18) does not account for the shearing resulting from the tire rotation that induces the lateral slip. Note that in this work the lateral slip is varied by the steering angle applied to the tire.

Recognizing this shortcoming in the state-of-the-art fast terramechanics model, a new surrogate model for SCM is developed in the next section.

II-C New Surrogate Model

This section presents the new terramechanics model developed as a fast surrogate for the SCM.

The new surrogate is obtained by replacing (17) and (19) with the following expressions, respectively.

[TABLE]

[TABLE]

with

[TABLE]

where $\Delta\delta$ is the step change in the steering angle. $k_{y}$, a parameter originally describing the shear displacement required to generate peak shear stress, is now empirically estimated as a function of the wheel states, i.e.,

[TABLE]

The lateral force acting on the vehicle is then determined as in (16). It should be noted that an additional term in (16) is often given representing the bulldozing force; however, simulations suggest this contribution is minimal for this application. Additionally, integrating the original nonlinear functions over the contact patch is a computationally demanding task. Therefore, the quadratic approximation proposed in [12] is adopted in the surrogate model. Furthermore, the modifications shown represent the lateral force acting on the vehicle frame, not the lateral forces in the tire reference frame.

Simulations covering the operating range of a notional military AGV are run to develop the modifying functions $g_{1}(\cdot)-g_{3}(\cdot)$. For each slip range of the clay simulation, as described in the Appendix, the simulations are run at 4 equispaced wheel loads, 5 slips, 5 translational velocities, and 7 sinkage exponents. Other terrain parameters are set to their nominal values, because only the sinkage exponent is selected as the parameter to be estimated due to the higher sensitivity of tire forces to sinkage exponent than other parameters [5, 19]. Table IV shows the range of inputs covered in the simulations. In these simulations, the inputs are held at constant values and the lateral slip is varied with a sinusoidal input. Following this, correction factors are determined for $\tau_{y}$, $j_{y}(\theta)$, and $k_{y}$ to match the output of (16) with each of the SCM simulations. Least squares curve fitting is then used to derive the relationship between the correction factors and the simulation inputs of Table IV, resulting in the modification functions $g_{1}(\cdot)-g_{3}(\cdot)$. To ensure the model was not subject to overfitting, over 1,500 independent validation simulations were performed as described in the following paragraphs.

Several parameters in the surrogate terramechanics model have distinct effects on the lateral force prediction and can be modified to achieve better agreement with SCM. For illustration purposes, the effects of each input on a sand terrain are shown in Fig. 3-5. The slope of the linear region can be set by modifying $k_{y}$ with $g_{3}(\cdot)$, the distance between the two curves can be set by adjusting $j_{y}$ with $g_{2}(\cdot)$, and the overall magnitude of the force can be adjusted with $g_{1}(\cdot)$. The effect of the wheel states on $g_{1}(\cdot)-g_{3}(\cdot)$ are as follows. The effect of wheel load can be captured with linear functions for $g_{1}(\cdot)-g_{3}(\cdot)$. Increased wheel load tends to increase the magnitude of the lateral force, increase the slope of the linear region, and increase the seperation between the top and bottom curves as shown in Fig. 3. The effect of longitudinal slip can be modeled by polynomials for $g_{1}(\cdot)-g_{3}(\cdot)$. As seen in Fig. 4, for positive slips, lower magnitude longitudinal slips tend to increase the slope of the linear region, while also causing a larger spread between the top and bottom curve. The effect of translational velocity can be captured with a power function for $g_{2}(\cdot)$ alone, because it has minimal effect on $g_{1}(\cdot)$ and $g_{3}(\cdot)$. Hence, $g_{1}(\cdot)$ and $g_{3}(\cdot)$ do not depend on translational velocity. As seen in Fig. 5, increased translational velocity tends to have little effect on the slope of the linear region, but reduces the hysteresis. Example formulations of $g_{1}(\cdot)-g_{3}(\cdot)$ for a clay terrain are given in the Appendix.

Once $g_{1}(\cdot)-g_{3}(\cdot)$ are determined, over 1,500 independent validation simulations are ran. The results of the surrogate model are shown in Fig. 6 (blue). Much better agreement is observed between the surrogate model and the SCM simulation compared to Fig. 2. It should also be noted that the surrogate model runs in 200-400 $\mu$s, which is an order of magnitude more efficient than what is reported for SCM and more suitable for online terrain estimation.

III Vehicle Models

Two vehicle models are employed in this work; a 11 DoF model to represent the plant and a 3 DoF bicycle model to be used as part of the estimator. This section describes these models.

III-A Plant Model

To represent the actual vehicle in the simulation-based validation of the proposed surrogate model and terrain estimator, the Chrono software is utilized to simulate the dynamics of a notional military vehicle as well as to implement the SCM [18]. For the purposes of this work, the vehicle is modeled with a double wishbone suspension, rack-pinion steering, 4 wheel drive, and a simple powertrain without a torque converter or transmission. This results in a 11 DoF vehicle model that is then combined with the SCM as the tire-terrain interaction model. The data received from the plant is then corrupted with Gaussian noise and serves as the measurement $y_{k}$ in (28) in Sec. IV. Table V lists the standard deviations used in the noise model for each state. Actual sensors typically offer lower noise levels; hence the chosen standard deviations represent a worse-case scenario to test the ability of the estimator [20].

III-B Bicycle Model

As part of the terrain estimation process that is detailed in Sec. IV, a vehicle model is needed to predict future vehicle states based on the tire forces from the surrogate model. For this work, a 3 DoF bicycle model with forward Euler integration is adopted, as it provides a proper level of fidelity while maintaining enough simplicity for short-horizon predictions [21]. The bicycle model takes on the following form:

[TABLE]

where the state vector, $z_{b}$, is defined as

[TABLE]

with $M_{t}$ being the vehicle mass, $I_{zz}$ being the vehicle’s yaw moment of inertia, and $L_{\text{f}}$ and $L_{\text{r}}$ being the distance from the vehicle’s center of gravity to the front and rear axles, respectively. Finally, $F_{\text{yf}}$ and $F_{\text{yr}}$ are the lateral forces generated from the front and rear tires acting on the vehicle body, as obtained from the terramechanics model.

IV Terrain Estimation

The terrain parameter to be estimated is chosen as the sinkage exponent $n$, because it has been shown to be the dominant parameter [5]. All other terrain parameters are assumed to be some nominal values based on the specific terrain type, which can be determined from terrain classification algorithms such as the ones described in [9, 22].

To estimate the unknown terrain parameter $n$, it is appended to the 3 DoF bicycle model in (24) with trivial dynamics. Here $n$ is given as a 2x1 vector to account for the front and rear tires. This is to mitigate the influence of unmodeled multipass effects and in the case of a discrete terrain change where the front tire and rear tire may operate on different terrains. The augmented state vector and state dynamics are given as

[TABLE]

Given the measurements of the vehicle states in (25), the augmented dynamics are utilized in an unscented Kalman filter (UKF) to estimate the augmented state vector in (26) including the sinkage exponent. It is worth noting that many other algorithms are available in the literature for nonlinear parameter estimation, including, but not limited to, extended Kalman filters, transitional Markov Chain Monte Carlo algorithms, and particle filters. Among these options the UKF is preferred in this work, because preliminary explorations suggest that the UKF offers a good balance between accuracy and computational speed for this application.

The UKF is composed of two general steps; a time update step and a measurement update step. Assume that a system is given in discrete time as:

[TABLE]

[TABLE]

where $z$ is the state, $y$ is the observation, and $v$ and $n$ are the process and observation noise, respectively. The functions $F(\cdot)$ and $H(\cdot)$ are nonlinear functions describing the dynamics and outputs. In this application, $z$ takes the form of the state vector in (26) and $F(\cdot)$ is obtained by discretizing the state equation in (26) using the forward Euler method. $H(\cdot)$ is given as the state vector in (25).

First, a set of $2L+1$ sigma points are created to capture the statistical distribution of the states, where $L$ is the dimension of the state vector $z$. The sigma points are determined as follows:

[TABLE]

where $\hat{z}$ is the mean value of $z$. $\lambda$ is a scaling parameter given as:

[TABLE]

where $\alpha$ is a tunable scaling parameter that typically takes a value between 0 and 1, $\kappa$ is another scaling parameter that is typically set to 0. At the time update step, the sigma points are propagated through the original nonlinear system as:

[TABLE]

The following weights are then calculated

[TABLE]

[TABLE]

[TABLE]

where $\zeta$ is set to 2 for Gaussian distributions. The statistics of the time update step are then given by:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, the measurement update step is given by the following set of equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The above process is a summary of the algorithm given in [23]. Intuitively, the process works by merging model-based predictions of the states with their measurements from the plant by exploiting the uncertainties associated with each to determine the best estimates of the states. For further discussion of UKF and details of its implementation, the reader is referred to [23, 24]. It also worth noting that while this work only focuses on estimating the dominant parameter $n$, other terrain parameters could be estimated simultaneously, as well. However, this would incur additional computational costs as the state space dimension increases, thus increasing the number of sigma points necessary in the UKF.

V Results and Discussion

In this section the performance of the terrain estimator is evaluated. The performance is evaluated from two different points of view: (1) the accuracy of the estimated sinkage exponent $n$, and (2) the accuracy of the predicted state trajectories of the vehicle . The former assesses the algorithm’s ability to find the true sinkage exponent, whereas the latter assesses the utility of estimating the sinkage exponent in the larger picture of predicting the future states of the vehicle. Note that if the assumed nominal values for the non-estimated parameters are not representative of the true terrain type, then the estimator may not necessarily converge to the true terrain parameter, as it will attempt to find a value that compensates for the errors in the non-estimated values and achieves the best prediction capability of the vehicle model. For the ultimate aim of more accurately predicting the future mobility capabilities of the vehicle, the second evaluation criterion is the more relevant one.

Simulations are performed utilizing Chrono’s SCM deformable terrain and the developed AGV model in Sec. III-A. Two terrains are considered including sandy loam and clay. Relevant terrain parameters are given in Table VI. The simulation subjects the AGV to sinusoidal steering commands, steering fully to the left and right over a three second period. The throttle is also varied with a sinusoidal command such that varying speeds are achieved. No braking command is given. The applied steering and speed profiles for the clay simulation are shown in Fig. 7. Throttle and steering commands of the same frequency are given in the sandy loam simulation, as well. Two remarks are in order. First, as seen in Fig. 7, no requirement on constraining the vehicle to low speeds (on the order of 10 cm/s) is made, which is in contrast to previous efforts [11]. This enables enhanced mobility, which is critical for military applications. Second, a sinusoidal steering input is selected to induce lateral dynamics for the vehicle. Since the bicycle model only utilizes the lateral forces acting on the vehicle, it is critical for the estimation that the vehicle operates in such a way that lateral dynamics are induced. Otherwise, the lack of information on the lateral dynamics leads to parameter variations having negligible effects on the output of the bicycle model. In other words, if $F_{\text{yf}}$ and $F_{\text{yr}}$ are zero, it is not possible to estimate terrain parameters based on lateral forces.

In all simulations, the simulation time step is set to 2 ms in Chrono. The purpose of these simulations is to determine the estimation algorithm’s accuracy and utility under different terrain conditions. Once the simulations are complete, noise is added to the outputs to simulate sensors, as discussed in Sec. III-A. The estimator is then run at a 12 ms time step and the simulated measurements are received at every 24 ms. The terrain parameter $n$ is initialized with a value off of the true terrain parameter used in the plant simulation. The remaining terrain parameters are set to their true values given in Table VI. Note that the true values are used here only to assess how closely the algorithm can converge to the true sinkage exponent.

Table VII displays the initial guess of the value of the sinkage exponent $n$, its converged estimate by the algorithm, and the error associated with the estimated parameter. The initial terrain parameter for sandy loam is chosen to be representative of Buchele (Michigan) sandy loam and the initial terrain parameter for clay is chosen to be representative of Thailand clay [25]. On both terrains, the percent error in the estimated terrain parameter is less than 4%, where the estimated value is taken to be the final value by the end of the simulation. Fig. 8 shows the estimated terrain parameters for the two considered terrains as time evolves. The differences between the converged and true terrain values can be due to model discrepancies between the high fidelity Chrono model and the 3 DoF bicycle model along with discrepancies arising from the reduced order terramechanics model. Nevertheless, the estimator converges within 10% of the estimated parameter within 5 seconds for both cases.

The peak computation time of the estimator is 10.5 ms and 7.5 ms for clay and sandy loam simulations respectively, thus demonstrating the potential to achieve real-time estimation. The platform running this estimation consists of 16 GB Memory and a single core 3 GHz Intel Core i7 processor.

While estimating the terrain properties accurately is a worthy goal in and of itself, it is more of interest to evaluate to what extent the estimations can improve the predictive capability of the bicycle model as motivated above. To accomplish this second evaluation, the bicycle model, with the terramechanics model parameterized by either the initial guess or the converged terrain parameter, is used to predict the vehicle states approximately 0.5, 2.5, and 5.0 seconds into the future for the clay case. After this time the vehicle states are reset to the true values received from Chrono. As such, this procedure mimics the operational procedure of a model predictive control approach, where a receding finite horizon optimal control problem is solved periodically with updated information available from sensors [3].

Table VIII depicts the mean squared errors (MSE) of the state estimates given by the 3 DoF bicycle model for both the case when the initial guess for the sinkage exponent for clay is used and the case when the converged estimate is used over the entire 32.89 s simulation. The model parameterized by the estimated terrain property yields significantly better predictions, especially at larger time horizons with order of magnitude reductions in MSE. Fig. 9 shows a portion of the simulation, using the $\sim$2.5s time horizon and depicting the true vehicle positions from the plant (blue solid line), and the predicted positions using the bicycle model with the initial guess of the sinkage exponent (black dotted line) and with the converged sinkage exponent (red dashed line). As can be seen, the converged value yields much more accurate predictions, thus demonstrating the ability of the estimator to significantly improve prediction fidelity. Similar results are also observed for the other terrains, but they are not reported here due to space limitations. This improvement in turn could lead to better performance in model predictive controllers, which is subject to future research.

While the results of the estimator are promising, there are two limiting assumptions of the proposed scheme. The first is that the estimator assumes the terrain is homogeneous. This is common among the approaches reported in the literature [11]. However, in reality, terrain parameters may be changing and evaluating the performance of the estimator in this scenario is subject to future work. The second limiting assumption is that SCM is treated as the ground truth and the surrogate terramechanics model is parameterized accordingly. The reality may be different than SCM. However, in that case, a similar procedure in developing the surrogate model could potentially be used by replacing the SCM simulations with an experimental single wheel test bed. Experimental validation of the developed terramechanics model is subject to future work.

VI Conclusion

This paper considers AGVs operating on deformable terrains with unknown terrain properties and develops a novel terrain estimation framework towards increasing the terrain-awareness of the AGV. In particular, the novelty of the framework is the development of a new surrogate terramechanics model for SCM and its use in conjunction with a bicycle model in a UKF. The results suggest that this new framework can estimate the dominant terrain parameter, namely the sinkage exponent, with high accuracy and high computational efficiency. It is therefore concluded that the framework is an important step towards achieving a good balance between estimation accuracy and computational speed. The results also show that the increase in the accuracy of the terrain parameter due to the developed estimation framework leads to a significant increase in the predictive accuracy of the bicycle model, especially for longer prediction time horizons. It is therefore concluded that the proposed framework could be useful to increase the performance of AGVs when they are controlled with model predictive schemes.

Future work includes evaluating the estimator on varying terrain conditions. The need for and ability of estimating multiple terrain parameters also needs to be investigated. It is also of interest to perform experimental validation of the estimator on an actual AGV, and investigate the utility of the estimator in a model predictive control framework for terrain-aware autonomous navigation.

Appendix

The following depicts the formulas for $g_{1}-g_{3}$ for a clay terrain. In this work the equations for $g_{1}-g_{3}$ are determined for four separate slip ranges for better agreement. The slip ranges are:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The dependencies for $g_{1}-g_{3}$ on each input is determined through curve fitting to simulation data as discussed in Sec. II-C. The below equations are valid for the slip range of $0\leq s<0.16$. For the lower curve of Fig. 6, the equations are given as:

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The same process yields the following equations for the upper curve of Fig. 6

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The equations for other slip ranges can easily be determined by repeating the curve fitting process on data in those ranges. Furthermore, it should be noted that these particular equations are only valid for the specific tire under consideration. Should these be used for a different tire, for example of different radius, the equations are no longer valid and the process would need to be repeated.