Abstract
This paper considers wheeled vehicles traversing unknown terrain, and proposes an approach for identifying the unknown soil parameters required for vehicle driving force prediction (drawbar pull prediction). The predicted drawbar pull can potentially be employed for traversability prediction, traction control, and trajectory following which, in turn, improve overall performance of off-road wheeled vehicles. The proposed algorithm uses an approximated form of the wheel-terrain interaction model and the Generalized Newton Raphson method to identify terrain parameters in real-time. With few measurements of wheel slip, i, vehicle sinkage, z, and drawbar pull, DP, samples, the algorithm is capable of identifying all the soil parameters required to predict vehicle driving forces over an entire range of wheel slip. The algorithm is validated with experimental data from a wheel-terrain interaction test rig. The identified soil parameters are used to predict the drawbar pull with good accuracy. The technique presented in this paper can be applied to any vehicle with rigid wheels or deformable wheels with relatively high inflation pressure, to predict driving forces in unknown environments.
Introduction
Wheeled vehicles traversing rough unknown terrain has many potential applications, including defence, agriculture, mining, and space exploration. Most wheeled vehicles still operate under driver supervision while carrying out tasks. Knowledge of the terrain characteristics helps the driver to have a better control of the vehicle.
From wheel-terrain interaction dynamics, it is seen that soil parameters play a vital role in determining vehicle drawbar pull which can, in turn, be utilized for developing traversability prediction criteria and traction control algorithms.
Research on wheel-terrain and track-terrain interactions has progressed since Bekker started pioneering this area (Bekker G., 1956), (Bekker G., 1969). Iagnemma and Dubowsky focused on the characteristics of a rigid wheel in deformable terrain for planetary rover applications (Iagnemma K. & Dubowsky S., 2002), (Iagnemma & et.al., 2004). In (Hutangkabodee S. & et.al., 2007), soil parameter identification for a tracked vehicle traversing unknown terrain is presented. Soil parameter estimation for a wheeled vehicle traversing deformable terrain was first carried out by Iagnemma and Dubowsky (Iagnemma K. & Dubowsky S., 2002). A Linear Least Square estimator was employed as an on-line identification method to identify two key soil parameters, cohesion (c) and internal friction angle (φ), using on-board rover sensors.
In this paper, for the first time, a method for identifying all five soil parameters is proposed for wheel-terrain interactions. The method is an extension from (Hutangkabodee S. & et.al, 2006) where three soil parameters were identified.
The Composite Simpson's Rule (CSR) is used to approximate the integrals of the original wheel-terrain interaction model, as they cannot be integrated analytically. By doing this, the speed of the algorithm increases by a factor of 9. The Generalized Newton Raphson (GNR) method is applied to the modified nonlinear wheel-terrain interaction model for soil parameter identification. This method has been shown to identify unknown parameters with high accuracy and rapid convergence (Zweiri Y.H. & et.al., 2004).
The experimental data from the wheel-terrain interaction test rig (Section 4) and the actual parameters of the experimental soils (obtained by methods described in Section 5) are used to validate the proposed algorithm. The test results show good identification accuracy for all the soil parameters except the pressure-sinkage parameters (S and n) for Garside 60 soil. The reason for the high errors is described in Section 6. The identified soil parameters are used to predict DP over the entire wheel slip range showing good agreement with the measured values. The predicted DP can be used for vehicle performance optimization, traversability prediction, traction control, and trajectory planning.
Governing Model
The analytical model of a wheeled vehicle interacting with an unknown deformable terrain is presented in this section, Fig. 1. In this figure, the wheel of a wheeled vehicle is considered rigid. For constant vehicle forward velocity, considering the force balance in the horizontal direction gives (Iagnemma K. & Dubowsky S., 2002):
where, DP is the drawbar pull applied to the wheel from the suspension system, θ is the wheel-terrain contact angle measured from vertical axis, θ1 is the angle at which the wheel first makes contact with the terrain, θ2 is the angle at which the wheel loses the contact from the terrain, σ is the radial stress normal to the wheel-terrain contact, τ is the shear stress tangent to the wheel-terrain contact, r is the wheel radius, b is the wheel tread, z is the sinkage of the wheel under the terrain surface, θ m is the angle at which the maximum radial stress occurs, σ1 is the radial stress profile between θ1 and θ m , and σ2 is the radial stress profile between θ m and θ2.

Free body diagram of a rigid wheel on deformable terrain (Iagnemma K. & et.al., 2004)
θ2 is assumed to be zero and the location of the point of maximum radial stress, θ
m
, is assumed to be in the middle between θ1 and θ2 (Iagnemma K. & et.al, 2004); therefore, θ
m
is equal to θ1/2. The model for soil parameter identification, (1), then becomes (Iagnemma K. & et.al, 2004):
where,
τ1(θ) and τ2(θ) are derived from
at σ(θ) = σ1(θ) and σ2(θ) respectively.
In (3), (4), and (5), k c , kφ, and n are the pressure-sinkage coefficients of the soil. These parameters have been derived for terramechanic applications and they characterize normal pressure from a vehicle with sinkage in the soil. c and φ are the soil cohesion, and internal friction angle, respectively. These two parameters characterize the strength of a soil when exposed to axial and shearing loads. K is the shear deformation modulus of the soil and characterizes displacement behavior of the soil under shear forces from a driven wheel.
All four integral terms of the DP model, (2), cannot be integrated analytically. The Composite Simpson's Rule (CSR) is applied to find an approximation to these integrals in order to facilitate the implementation of the identification algorithm and increase the identification speed.
For the integral
The drawbar pull (DP) based on the Composite Simpson's Rule (CSR) approximation and the DP based on numerical integration are shown in Fig. 2. It is observed from this figure that the CSR gives a very accurate approximation of the DP model with an rms error = 0.22 %. In deriving the DP curves in Fig. 2, it is found that while the numerical integration takes 0.045 s, the CSR approximated model takes only 0.005 s. Therefore, this new approximate DP model is more suitable for soil parameter identification in real-time.

Comparison of the DP using the CSR approximation and numerical integration
The pressure-sinkage coefficient, k c /b + kφ, is a constant, for constant b values. Since in current applications the wheel tread b of a vehicle is constant, a lumped pressure-sinkage coefficient, S = k c /b + kφ, is used in (2) and will be treated as a single soil parameter.
The proposed soil parameter identification algorithm for wheeled vehicles traversing unknown terrain is based on the Generalized Newton Raphson (GNR) method. This method is well suited for real-time parameter identification since it has been shown to identify unknown parameters with high accuracy and rapid convergence (Zweiri Y.H. & et.al., 2004).
The concept of the GNR method is based on the Newton Raphson (NR) method. The NR method works by iteratively modifying an initial estimate and, after a number of iterations, arrives at a converged solution. Commonly, the iterative process is halted when the error between the identified value and the current value falls below a predefined threshold. The main difference between the GNR and the NR methods is that the NR method employs the same number of equations as the number of unknowns whereas the GNR method uses more equations than the number of unknowns.
The GNR method has the same advantages as those of the NR method, namely its ability to identify unknown individual parameters, its robustness to noise and its fast speed of convergence. In addition, the robustness of the GNR method is even better than that of the NR method (Zweiri Y.H. & et.al., 2004) because it considers more measured data points than the number of unknowns. The GNR method also allows us to use all the measured data available (instead of having to pick 5 data points from maybe 10 available data points as in the NR method case), and thus helps further improve the identification accuracy. However, it is noted that if the initial estimate is selected at a point where the derivative of the function is relatively small, the convergence speed may be slow (Faires J.D. & Burden R, 2003).
The GNR method is implemented to identify the soil parameters for vehicle wheel-terrain interaction dynamics, Fig. 3. Here, vector p comprises five parameters, the cohesion (c), the internal friction angle (φ, the shear deformation modulus (K), the lumped pressure-sinkage coefficient (S), and the pressure-sinkage exponential (n). The vector x consists of three measured signals, drawbar pull (DP), sinkage (z), and wheel slip (i).

Diagram showing implementation of the GNR method for soil parameter identification
The Composite Simpson's Rule modified wheel-terrain interaction model of (2), and the measurement vector x are used to identify the unknown soil parameters. To find p, q independent equations are required (q ≥ 5). These equations are generated by evaluating the function f at q different time samples (t1, t2, …, t
q
). This results in q nonlinear equations expressed in matrix form as:
where, p = [c, φ, K, S, n] T , and x = [DP, z, i] T .
Let
where,
The wheel-terrain interaction test rig used in this study is shown in Figs. 4 and 5.

Wheel-terrain interaction test rig (1): where, A is the wheel assembly, B is the motor-gearbox-encoder assembly for driving carriage, C is the chain, D are the shafts connected from the top of wheel assembly, E is the carriage, F is the potentiometer, G is the test rig frame, and H is the experimental soil

Wheel-terrain interaction test rig (2): where, I is the rigid wheel, J is the wheel motor-gearbox-encoder, K is the Gamma 6-Axis force/torque transducer, L is the moving part of the potentiometer for measuring distance (called “wiper”), and M is the connecting rod for potentiometer's wiper and wheel assembly.
This test rig is developed to test the algorithm proposed in Section 3 in a controlled environment allowing us to acquire all needed parameters in a simple and straightforward manner, and to test the feasibility and effectiveness of the soil parameter identification algorithm. In particular, it is difficult to measure slip parameters on a real vehicle manoeuvring through terrain. Based on current developments in slip estimation (Le A.T., 1999), (Song Z. & et.al., 2005), (Ohno K. & et.al., 2003), (Barshan B. & Durrant-Whyte H.F., 1995), (Chhaniyara S. & et.al., 2006), we assume slip to be measurable on-line enabling the implementation of the proposed soil parameter identification process. With this test rig, the wheel slip, i, the wheel sinkage, z, and the drawbar pull, DP, can be measured very accurately. This allows the soil parameter identification algorithm to be validated in an effective way.
In practice, the signals needed for soil parameter identification in real wheeled vehicles can be acquired as follows. To measure the drawbar pull (DP) of the vehicle, a single-axis force sensor can be installed between the wheel and the vehicle suspension to measure the resultant force in the vehicle heading direction. The vehicle sinkage, z, can be acquired by measuring the gap between axle of driven wheel and ground. Ultrasonic sensors (facing downwards) attached to the vehicle suspension system can be employed for this purpose.
Slip is estimated by employing a sliding-mode observer whose inputs are wheel angular speeds and vehicle heading speed (Song Z. & et.al., 2005). The wheel angular speeds can be measured by attaching either shaft encoders or tachometers onto vehicle wheels' shafts. The vehicle heading speed can be obtained using several techniques such as Global Positioning System (GPS), Inertial Navigation System (INS) (Barshan B. & Durrant-Whyte H.F., 1995), and Optical Flow (Chhaniyara S. & et.al., 2006).
The test rig consists of a driven wheel mounted on a horizontal driven carriage. The wheel is not vertically constrained and can move freely in the vertical direction. The wheel assembly is connected to the carriage assembly using 2 shafts. During experimentation, two motors - one driving the wheel and another driving the carriage - are operated simultaneously at different speeds to generate wheel slip. By varying the differential speeds between the wheel and the carriage motors, various controlled wheel slips can be generated. The test rig is equipped with various sensors to acquire all the measured signals needed for soil parameter identification. Two optical encoders are used to measure the angular speeds of the wheel motor and chain motor. The wheel slip can then be calculated using (9).
where, i is the wheel slip, r w is the wheel radius, r c is the chain sprocket gear radius, ω w is the wheel angular velocity, ω c is the chain sprocket gear angular velocity.
A Gamma ATI 6-Axis force/torque (F/T) transducer mounted between the wheel assembly and the horizontal carriage is used to acquire the drawbar pull, DP. The sinkage (z) of the wheel is measured by a potentiometer. The wheel used in the experiment is 0.1 m in width and 0.2 m in diameter. The test rig is fixed to a 1 m × 1 m × 0.5 m frame. The experimental soil is contained in a 0.3 m × 0.85 m × 0.1 m soil box. Experiments were performed on two types of soils - Garside 60 soil, and Garside iron sand. These soils were purchased from the Garside Sands Company, UK.
Parameters c and φ
To obtain the cohesion, c, and internal friction angle, φ, of the experimental soils, the Direct Shear Test is conducted (Terzaghi K., 1948). The test is performed on three specimens from a relatively undisturbed soil sample. A specimen is placed in a shear box which possesses two stacked rings to hold the sample. A normal load is applied vertically to the specimen (creating a certain normal stress to the specimen). The upper ring is pulled laterally until the sample fails and the maximum shear stress is then recorded.
Another two different normal loads are applied to another two specimens for the same soil type. The specimens are again sheared until failed and the corresponding maximum shear stresses are stored. The results for these three specimens (normal stress-shear stress values) are plotted on a graph with shear stress as the y-axis and normal stress as the x-axis, Fig. 6. The y-intercept of the experimental curve represents c, and the slope of the curve is φ of the tested soil.

Shear stress vs. normal stress plot and definitions of c and φ
The pressure-sinkage parameters k
c
, kφ, and n of the experimental soils are described in the following formula:
where, p is the pressure exerted to the soil, b is the diameter of a circular disc used as the contact area to the soil (or the width of a rectangular plate), z is the corresponding sinkage of the disc into the soil.
Equation (9) characterizes the sinking behavior of a soil when vertically pressed with a circular disc or a rectangular plate. To obtain all three parameters (k
c
, kφ, and n), two discs/plates with different diameter/width are required. Equations (10) and (11) show two sets of equations employed for this purpose (with logarithmic operation applied on both sides of the equation).
where, b1 and b2 are the diameters of circular discs 1 and 2 used in this experiment, respectively.
These two equations can be shown schematically in Fig. 7. To obtain these two curves, two different sizes of circular discs are pressed to sink dynamically into the soil whose pressure-sinkage cofficients (k c , kφ, and n) are to be measured. The reason the pressure and its corresponding sinkage have to be measured in the dynamic manner is that it simulates real circumstances where a vehicle is moving on a soil (not like in civil engineering case where static load is placed on a certain soil and corresponding sinkage is measured when an object stops to sink). k c , kφ, and n are not valid in civil engineering field. They are terms created for soil in terramechanic applications.

Pressure-sinkage curves for Equations (10) and (11)
The traditional way to carry out the pressure-sinkage experiment is to use a hydraulic system with pressure gauge installed (Bekker G., 1956), (Wong J.Y., 2001). However, due to high cost associated with the system, the Screw Press system is invented to acquire parameters k c , kφ, and n, Fig. 8.

Screw Press system developed for pressure-sinkage experiment, where, A is the screw press, B is the screw press's handle, C is the potentiometer fixed with the static part of the screw press, D is the experimental soil, E is the S-shape strain gauge, and F is the circular disc.
With the soil normal pressure (exerted from the disc) measured by the S-shape strain gauge (label E, Fig. 8) and the corresponding disc's sinkage measured by the potentiometer or linear position transducer (label C, Fig. 8), parameters k c , kφ, and n can be obtained by solving Equations (10) and (11).
Two methods described in (Wills B. M. D., 1963) are employed for measuring K. The first method uses annular torsional shear test (Hvorslev M. J., 1939). The apparatus for this test is driven by a hydraulic ram on an experimental soil. The values of applied torque and corresponding soil deformation are recorded using a special indicator and they are used for K calculation (Wills B. M. D., 1962). The second method involves using translational rigid track. The track is pulled by a hydraulic ram, and the shear force and corresponding soil deformation are measured from the oil pressure and movement of the piston, respectively (Wills B. M. D., 1962).
Experimental Results
In this section, the soil parameter identification results are presented based on experimental data taken from the test rig described in Section 4. From all the measured data obtained from experiments, the entire sets of measured DP, z, and i are used for soil parameter identification based on the algorithm developed in Section 3. There are 20 and 17 sets of data for Garside 60 soil and Garside iron sand respectively, acquired from wheel-terrain interaction test rig. The measured data for both soils were curve-fitted to a polynomial spline before being used as inputs to the algorithm. Tables 1 and 2 show the identification results for all 5 soil parameters for Garside 60 soil and Garside iron sand, respectively.
Identification results of 5 soil parameters for Garside 60 soil
Identification results of 5 soil parameters for Garside 60 soil
Identification results of 5 soil parameters for Garside iron sand
The identification algorithms were run on an Intel Pentium(R) 4 processor with a 2.80 GHz CPU and 1.00 GByte of RAM.
K varies from 0.01 m (firm sand) to 0.025 m (loose sand) for sand-like soils (Wong J.Y., 2001). The parameter K of both experimental soils falls within this range, as they are sands with varying degrees of compaction (firmness). In this study, since K is very difficult to measure and requires complex measuring devices (as described in Section 5.3), its value is assumed to be between 0.01 m and 0.025 m. Values of c and φ for both soils were obtained by the Direct Shear Test (Section 5.1), whereas their k c , kφ, and n were obtained from pressure-sinkage experiments (Section 5.2).
From Table 2 (using test rig data for Garside iron sand), the identification accuracies are relatively high with identification errors less than 12.52 % for all 5 parameters. For parameter K, ‘In range’ shown in percentage error column means that the identified K is within range of the actual K (0.01 m – 0.025 m). From Table 1 (for Garside 60 soil), the identification errors are quite high for parameters S and n (54.72 % and 25.35 %, respectively); however, they can be used to predict the drawbar pull with good accuracy (7.25 % rms error), Fig. 7 and Table 1. ‘N/A’ shown in percentage error column (Table 1) means that error for c cannot be computed since its actual value is equal to 0. However, its identified value is considered close to its actual value (the practical range of c for real-world soils is 0 Pa ≤ c ≤ 68500 Pa (Wong J.Y., 2001)).
Figures 9 and 10 are used to explain the reason for high identification errors of S and n. As illustrated in these figures, with a given set of values for c, φ, and K (c = 1000 Pa, φ = 34 degree, and K = 0.015 m), using different sets of values for S and n gives similar DP plots. When compared to the DP calculated using S = 500000 N/mn+2 and n = 0.64, the rms error of the DP calculated using S = 50000 N/mn+2 and n = 0.2 is 0.84 % (Fig. 9), while that of the DP calculated using S = 8800000 N/mn+2 and n = 1.2 is 1.15 % (Fig. 10). It can, therefore, be concluded that a slight deviation of DP inputs to the algorithm would yield different set of results for S and n. In reality, due to the random nature of the measurement noise, the measured DP will result in inaccuracies of identified pressure-sinkage parameters, S and n.

Plot of DP using c = 1000 Pa, φ = 34 degree, K = 0.015 m with different pairs of S and n (case 1 when using low-valued n)

Plot of DP using c = 1000 Pa, φ = 34 degree, K = 0.015 m with different pairs of S and n (case 2 when using high-valued n)
The identification times are 0.516 s and 0.344 s for Garside iron sand and Garside 60 soil, respectively. The execution time can be further reduced if the code is optimized. In this paper, the code is created using Matlab 6.5 running in command mode; the execution time can be significantly reduced, if an optimized code is executed. Therefore, the proposed algorithm has promising potential for on-line soil parameter identification for off-road wheeled vehicles.
The identification results from Tables 1 and 2 are used to predict the DP and compared with the experimental DP from the test rig, Figs. 11 and 12. From these figures, the rms errors are 7.25 % and 5.50 % for Garside 60 soil and Garside iron sand, respectively. This shows relatively good prediction accuracy of DP over the entire slip range using the identified soil parameters.

Comparison between the experimental DP and that predicted from the identified ‘c, φ, K, S, n’ for Garside 60 soil

Comparison between the experimental DP and that predicted from the identified ‘c, (, K, S, n’ for Garside iron sand
Next, a robustness test is carried out to examine whether the proposed algorithm will result in converged solutions when different initial conditions (p0) are used. Table 3 gives the ranges of initial conditions that allow successful convergence. It is noted that values in brackets in this table are the practical range of parameters for real-world soils (Wong J.Y., 2001).
The range of initial estimates of soil parameters that produces the converged solutions for both experimental soils
The algorithm shows good robustness to initial c and φ since any value chosen from their practical range results in converged solutions. However, for the other three parameters, it is recommended that initial K be chosen towards its practical lower bound (0.006 m ≤ K ≤ 0.014 m), initial S be chosen towards its practical upper bound (140000 N/mn+2 ≤ S ≤ 6321000 N/mn+2) and initial n be chosen towards its practical upper bound (0.15 ≤ n ≤ 1.6). It is noted that the practical lower and upper limits of parameter S varies upon the wheel's width of the wheeled vehicle employed.
To sum up, by measuring DP and z samples at 5 or more slip values, soil parameters can be identified and used to predict the vehicle DP over the entire range of slips. These identified soil parameters can also be sent to following vehicles to improve their navigation on terrains ahead. Hence, accurate DP prediction of a wheeled vehicle over a particular terrain can be achieved. Consequently, based on accurate prediction of vehicle DP, the proposed algorithm can potentially be advantageous for traversability prediction and traction control for off-road wheeled vehicles.
This paper investigates soil parameter identification for wheeled vehicles traversing unknown terrain based on a wheel-terrain interaction dynamic model and sensor feedback (drawbar pull, wheel slip, and wheel sinkage). The paper presents a new real-time approach for identifying all the unknown soil parameters required for drawbar pull prediction based on the Generalized Newton Raphson (GNR) method. This is the first time a technique for identifying all the unknown soil parameters for wheeled vehicles is proposed. The experimental data from the wheel-terrain interaction test rig are used to validate the feasibility and effectiveness of the approach. The results show that, despite high identification errors for some parameters, the identified soil parameters can be used to predict the drawbar pull with good accuracy. The range of initial estimates for each parameter that give converged solutions is also proposed. With terrain parameters identified by the proposed algorithm, the vehicle drawbar pull for off-road wheeled vehicles can be predicted for all driving conditions. The identified terrain parameters can also be forwarded to vehicles following behind to improve their navigation on that terrain. The technique can be applied to any wheeled vehicle with rigid wheels or deformable wheels with relatively high inflation pressure.
Future work will focus on the outdoor testing with a real wheeled vehicle or wheeled mobile robot.
