Vibration-Based Recognition of Wheel–Terrain Interaction for Terramechanics Model Selection and Terrain Parameter Identification for Lugged-Wheel Planetary Rovers

Identifying terrain parameters is important for high-fidelity simulation and high-performance control of planetary rovers. The wheel–terrain interaction classes (WTICs) are usually different for rovers traversing various types of terrain. Every terramechanics model corresponds to its wheel–terrain interaction class (WTIC). Therefore, for terrain parameter identification of the terramechanics model when rovers traverse various terrains, terramechanics model switching corresponding to the WTIC needs to be solved. This paper proposes a speed-independent vibration-based method for WTIC recognition to switch the terramechanics model and then identify its terrain parameters. In order to switch terramechanics models, wheel–terrain interactions are divided into three classes. Three vibration models of wheels under three WTICs have been built and analyzed. Vibration features in the models are extracted and non-dimensionalized to be independent of wheel speed. A vibration-feature-based recognition method of the WTIC is proposed. Then, the terrain parameters of the terramechanics model corresponding to the recognized WTIC are identified. Experiment results obtained using a Planetary Rover Prototype show that the identification method of terrain parameters is effective for rovers traversing various terrains. The relative errors of estimated wheel–terrain interaction force with identified terrain parameters are less than 16%, 12%, and 9% for rovers traversing hard, gravel, and sandy terrain, respectively.


Introduction
Planetary rovers need to traverse rough terrains during planetary exploration missions.Wheel-terrain interaction has substantial influence on the mobility of rovers [1,2].The typical surfaces of Mars consist of hard, sandy, and gravel surfaces [3].Rovers can move quickly on hard, flat, and non-slippery surfaces.However, they are typically required to traverse over sandy areas, where they suffer from wheel sinkage, or over gravel surfaces where violent vibrations may happen.Hazards originating from the terrain surface itself are called non-geometric hazards.Knowledge of the terrain class and other properties can help improve the mobility of the rover by adjusting its path or control strategy in advance.Recognition of terrain class, terramechanics modeling, and terrain parameter identification are essential to optimize the control of high-mobility rovers.Research progress and existing problems related to terrain classification and parameter identification are discussed in this section.
Sensors 2023, 23, 9752 2 of 32 Terrain classification can be divided into exteroceptive and proprioceptive terrain classification.Exteroceptive terrain classification mainly applies vision or range data [4][5][6][7][8].For example, Otsu et al. [4] used semi-supervised algorithms instead of supervised algorithms for visual classification systems.Support vector machine (SVM) is used to classify data from different sensors.For the excellent performance of deep learning in object classification and recognition, some studies have utilized deep learning algorithms for terrain classification of planetary rovers, achieving high classification accuracy [6,7].However, the vision-or laser-based classification methods are usually related to the topmost terrain surface, which may not be the load-bearing surface.An example of this is the Mars terrain, which is covered by a thin layer of drift sand.
Proprioceptive terrain classification has focused on vibration-based terrain classification methods [9][10][11][12][13][14]. Brooks and Iagnemma [9] extracted the power spectral density (PSD) of vibration acceleration to classify the terrain into gravel or sand.In a later work, vibration, vision, and wheel-terrain interaction force were combined to construct a self-supervised terrain classification framework for planetary rovers [10].However, the terramechanics model used in their paper is not suitable for calculating the interaction force of the wheel and hard/gravel terrain.The terrain parameters used in the terramechanics model are also unidentified.Vicente et al. [11] extracted eleven statistical features of vibration accompanied with the PSD of vibration acceleration to form a new feature vector.An extreme learning machine was used to classify terrains based on the feature vector.Some studies applied deep learning algorithms to classify planetary terrain based on proprioceptive sensors [12,13].Those methods can reduce classification accuracy caused by inappropriate feature extraction and selection.These methods neglect that the variation of the wheel speed (over a given terrain) changes the amplitude and frequency of the measured wheel vibration signal.If rovers move on the same class of terrain at different speeds, the class of terrain may be incorrectly classified.One solution is to increase the number of trained speeds to obtain vibration features at various speeds [14].However, there is a drawback that numerous samples at the rover's multiple moving speeds need to be obtained.Another solution is to extract speed-independent vibration features of the wheel that are proposed in this paper.
Terrain parameters are useful for planetary rovers to maximize the wheel traction or minimize the power consumption [15].Several methods have been reported to estimate terrain parameters [16][17][18][19].Iagnemma et al. [16] presented a linearized terramechanics model for estimating the internal cohesion and friction angle of soil.However, many other important terrain parameters such as the cohesive modulus, friction modulus, and sinkage exponent in terramechanics models have not been identified.Hutangkabodee et al. [17] identified the parameters, including the internal friction angle, shearing deformation modulus, and lumped pressure sinkage coefficient, using the Newton-Raphson method based on a simplified terramechanics model.Ding et al. presented an improved terramechanics model with consideration of the slip-sinkage and lug effects.Soil parameters were divided into three independent sets by decoupling a wheel-soil terramechanics model.This model can be used to identify Lunar terrain parameters with the Yutu-2 rover [18].However, these methods cannot be applicable to rovers traversing firm soil because the assumed terramechanics models do not cover interaction between the rover wheels and hard terrain (HT).The soil of HT does not experience a shear deformation, and there is no shear stress in the wheel-terrain interaction.The soil cohesion coefficient (c) or shearing deformation modulus (K) is not part of the shear characteristic parameters in the shear terramechanics model of the interaction between a wheel and HT.To address this problem, Yang et al. [19] set parameters K and c equal to zero to make a terramechanics model of the interaction between a wheel and soft terrain (ST) suitable for the interaction between the wheel and HT.It is necessary to carry out terrain classification to determine parameters such as K and c in a terramechanics model when calculating interaction force or torque.
Conventional terrain classification classifies terrain classes according to their components to help rovers and mobile robots adjust their control strategy and path planning Sensors 2023, 23, 9752 3 of 32 through which a qualitative judgment of the terrain properties is acceptable.However, an autonomous planetary rover should "know" the terrain properties to calculate the wheel-terrain interaction force to maximize the wheel traction or minimize the power consumption and improve the traversability of the rover during its operational phase.Despite efforts in this regard, the terrain parameter identification of various terrain classes and terramechanics model switching for the interaction between the wheel and soft or hard terrain are still issues to be researched for planetary rovers.
As mentioned above, the following problems are still associated with research works on terrain classification and parameter identification for planetary rovers: (1) a suitable terramechanics model needs to be switched to for rovers moving on various types of terrain.The recognition of wheel-terrain interaction class (WTIC) has been rarely studied.
(2) Parameter identification of terramechanics models in various types of terrain (not just sand) needs to be studied.And (3) reported vibration-based speed-independent terrain classification methods require numerous training samples at multiple speeds, which is not feasible for planetary rovers.To solve these problems, the following need to be investigated: (1) how to extract speed-independent vibration features and recognize WTIC to determine the corresponding terramechanics model, especially in the regions where the terramechanics model needs to be switched, and (2) how to identify parameters in the terramechanics model for a planetary rover traversing various types of terrain.
In this paper, the speed-independent vibration features are extracted to recognize WTIC first.Terrain parameters in the terramechanics model that is switched to corresponding to the recognized WTIC are then identified.This paper comprehensively utilizes various sensors installed on the planetary rovers, including accelerometers, cameras, and six-dimensional force/torque sensors, to achieve the recognition of WTIC and identification of terrain parameters.The accelerometer is used to collect wheel vibration acceleration data.Speed-independent vibration features are extracted from the acceleration data to recognize WTIC and select the corresponding terramechanics model.The camera can obtain wheelterrain interaction images to realize the estimation of wheel sinkage and slip ratio using machine vision technology.The wheel-terrain interaction force information of the rover wheel is acquired using the six-dimensional force/torque sensor.Then, the wheel-terrain interaction force, wheel sinkage and slip ratio are input into the terramechanics model to identify terrain parameters.
The main contributions of this paper are as follows: (1) corresponding to the three wheel-terrain interaction classes (WTICs), this paper proposes three vibration models of lugged wheels and simplifies them to be approximated as periodic functions, showing the difference in wheel vibration frequency and amplitude for different WTICs.Through the dimensionless processing of the amplitude, frequency, mean, and mean square value of the vibration acceleration, four speed-independent vibration features are extracted to reflect the influence of the terrain on the vibration of the wheel.(2) A vibration-based WTIC recognition method without requirement for training data is proposed.Results of single-wheel test experiments show that the simple and effective method has recognized 292 samples of 300 samples correctly, so the recognition accuracy is about 97%.(3) Corresponding to the recognized WTIC, the terramechanics model is switched by determining whether the parameters of the shearing deformation modulus (K) and soil cohesion coefficient © in the model are 0. Terrain parameters in the terramechanics model are identified using the least squares method for rovers traversing various types of terrain.Experiment results obtained using a Planetary Rover Prototype moving on three terrain classes show that the proposed vibration-based WTIC recognition method is effective and has recognized 411 samples of 432 samples correctly, and the terrain parameter identification for various terrains is also effective.
The remainder of this paper is organized as follows.Section 2 discusses the vibration models of the wheel for different WTICs.Based on these models, speed-independent features extracted from the vibration acceleration are described in Section 3. Section 4 presents a method for WTIC recognition by analyzing the range of vibration features for Sensors 2023, 23, 9752 4 of 32 different WTICs.In Section 5, a method is proposed for identifying parameters of the terramechanics model in various types of terrain.The method switches the corresponding terramechanics model according to the recognized WTIC and then uses the model to identify terrain parameters.Section 6 discusses the experiment results of the WTIC recognition and parameter identification for a Planetary Rover Prototype.Finally, some concluding remarks are provided in Section 7.

Vertical Dynamics Model of Wheel Motion
A lugged wheel can be modeled as a mass-spring-damper system shown in Figure 1 [20], which is composed of the mass m w that the wheels of the rover need to load, a spring with a stiffness coefficient k w , and a damper with a damping coefficient c w .Here, y(t) is the vertical output displacement, and q(t) is the base excitation, which is a time-varying vertical displacement of the lowest point of the wheel rim.In addition, Y(jω) and Q(jω) are the Fourier transforms of y(t) and q(t), respectively, where Y(jω) = H(jω)Q(jω).The frequency response function H(jω) is as follows: The remainder of this paper is organized as follows.Section 2 discusses the vibration models of the wheel for different WTICs.Based on these models, speed-independent features extracted from the vibration acceleration are described in Section 3. Section 4 presents a method for WTIC recognition by analyzing the range of vibration features for different WTICs.In Section 5, a method is proposed for identifying parameters of the terramechanics model in various types of terrain.The method switches the corresponding terramechanics model according to the recognized WTIC and then uses the model to identify terrain parameters.Section 6 discusses the experiment results of the WTIC recognition and parameter identification for a Planetary Rover Prototype.Finally, some concluding remarks are provided in Section 7.

Vertical Dynamics Model of Wheel Motion
A lugged wheel can be modeled as a mass-spring-damper system shown in Figure 1 [20], which is composed of the mass mw that the wheels of the rover need to load, a spring with a stiffness coefficient kw, and a damper with a damping coefficient cw.Here, y(t) is the vertical output displacement, and q(t) is the base excitation, which is a time-varying vertical displacement of the lowest point of the wheel rim.In addition, Y(jω) and Q(jω) are the Fourier transforms of y(t) and q(t), respectively, where Y(jω) = H(jω)Q(jω).The frequency response function H(jω) is as follows: ) The modulus of H(jω) can be expressed as: For a lugged wheel on the rover, the spring stiffness kw is very big.The mass mw is assumed as 29 kg, the same as the mass of Spirit rover equally distributed to each wheel.The rover always travels at a low speed, for example, the speed of the Spirit rover is less than 50 mm/s.As a result, the natural frequency ωn is much bigger than the frequency v f of base excitation that is an approximate periodic function for a lugged wheel, as shown The modulus of H(jω) can be expressed as: where ξ = c w /2 √ k w m w is the damping ratio, λ = f v /ω n is the frequency ratio, f v is the frequency of base excitation q(t), and ω n = √ k w /m w is the natural frequency of the undamped free mass-spring-damper system.According to [20], |H(jω)| is approximately equal to 1 when λ is less than 0.75.
For a lugged wheel on the rover, the spring stiffness k w is very big.The mass m w is assumed as 29 kg, the same as the mass of Spirit rover equally distributed to each wheel.The rover always travels at a low speed, for example, the speed of the Spirit rover is less than 50 mm/s.As a result, the natural frequency ω n is much bigger than the frequency f v of base excitation that is an approximate periodic function for a lugged wheel, as shown later in this paper.Therefore, λ is less than 0.75 and |H(jω)| is approximately equal to 1. Hence, According to the properties of a Fourier transform, the following approximation can be obtained.y(t) ≈ q(t).

of 32
Therefore, the vertical output displacement and base excitation of the rigid lugged wheel can be regarded as equal in value, i.e., y(t) = q(t).Thus, in the following section, the vibration base excitations of the wheel moving on various terrain classes are analyzed to obtain the wheel vibration output displacements.

Vibration Models for Wheels Moving on Different Terrains in Wheels' Center Plane
Wheels with lugs parallel to the wheel axis are usually used on planetary rovers, such as Spirit rover's wheels.This type of wheel is studied in this paper.The vibrations of this type of wheel result from lug motion and terrain roughness/hardness.Terrains are mainly divided into soft and hard terrain.Both of them can be divided into low-unevenness and high-unevenness terrain.When the planetary rover travels on different classes of terrain, the terramechanics model is different and needs to be switched by setting the values of terrain parameters K and c in the model.In order to determine the values of K and c, wheelterrain interactions are divided into three classes as shown in Figure 2: lugs entering into soil (LEIS), all lugs contacting terrain surface (ALCT), some lugs contacting terrain surface (SLCT).When the wheel moves on low-unevenness soft terrain and high-unevenness soft terrain, every lug can enter into the soil.When the wheel moves on low-unevenness hard terrain, every lug can contact the terrain surface.When the wheel moves on highunevenness hard terrain, the same lugs can contact the terrain surface while other lugs cannot contact with the terrain for a wheel rotating 360 • .The terrain is classified into three groups: the first class of terrain (T1), which corresponds to WTIC 1: ALCT, prevents each lug from entering the soil during wheel movement, although each lug can be in contact with the terrain.The second class of terrain (T2), which corresponds to WTIC 2: SLCT, makes it impossible for the lugs to enter the soil during the wheel movement, and some lugs cannot be in contact with the terrain.The third class of terrain (T3), which corresponds to WTIC 3: LEIS, refers to the terrain when the wheels transverse it, and where the lugs can enter the soil.The correspondence between terrain classes and the WTICs is shown in Table 1.
According to the properties of a Fourier transform, the following approximation can be obtained.
( ) ( ) Therefore, the vertical output displacement and base excitation of the rigid lugged wheel can be regarded as equal in value, i.e., ( ) ( ) y t q t = .Thus, in the following section, the vibration base excitations of the wheel moving on various terrain classes are analyzed to obtain the wheel vibration output displacements.

Vibration Models for Wheels Moving on Different Terrains in Wheels' Center Plane
Wheels with lugs parallel to the wheel axis are usually used on planetary rovers, such as Spirit rover's wheels.This type of wheel is studied in this paper.The vibrations of this type of wheel result from lug motion and terrain roughness/hardness.Terrains are mainly divided into soft and hard terrain.Both of them can be divided into low-unevenness and high-unevenness terrain.When the planetary rover travels on different classes of terrain, the terramechanics model is different and needs to be switched by setting the values of terrain parameters K and c in the model.In order to determine the values of K and c, wheel-terrain interactions are divided into three classes as shown in Figure 2: lugs entering into soil (LEIS), all lugs contacting terrain surface (ALCT), some lugs contacting terrain surface (SLCT).When the wheel moves on low-unevenness soft terrain and high-unevenness soft terrain, every lug can enter into the soil.When the wheel moves on low-unevenness hard terrain, every lug can contact the terrain surface.When the wheel moves on high-unevenness hard terrain, the same lugs can contact the terrain surface while other lugs cannot contact with the terrain for a wheel rotating 360°.The terrain is classified into three groups: the first class of terrain (T1), which corresponds to WTIC 1: ALCT, prevents each lug from entering the soil during wheel movement, although each lug can be in contact with the terrain.The second class of terrain (T2), which corresponds to WTIC 2: SLCT, makes it impossible for the lugs to enter the soil during the wheel movement, and some lugs cannot be in contact with the terrain.The third class of terrain (T3), which corresponds to WTIC 3: LEIS, refers to the terrain when the wheels transverse it, and where the lugs can enter the soil.The correspondence between terrain classes and the WTICs is shown in Table 1.The origin of coordinate system Σo-xy is located at the wheel center when the wheel is stationary.The x-axis is parallel to the horizontal direction, and the y-axis is opposite the gravity direction.Σo-xy does not move with the wheel.The angle velocity of the wheel is ω w , l is the lug height, n l is the number of lugs, and α nl (α nl = 2π/n l ) is the angle between two adjacent lugs.According to the number of lugs on the wheel, the wheel can be divided into two classes: olig-lugged and multi-lugged wheels.An olig-lugged wheel is defined as the lugged wheel whose rim can contact the terrain surface when traversing hard terrain (HT).A multi-lugged wheel, which is often used on rovers, is defined as the wheel equipped with a large number of lugs that keep its rim from making contact with the terrain surface when traversing HT.

Vibration Model of Wheels Moving on Flat Hard Terrain
Flat hard terrain is regarded as the reference terrain.The vertical vibration of a wheel traversing the reference terrain is defined as the vertical reference vibration.
Figure 3 shows the movement phases for a wheel rotating around a supporting lug on a reference terrain.Here, θ (−α nl /2 < θ ≤ α nl /2) is the angle from the direction of the support lug to the gravity direction at the initial time when the wheel starts to move in the anticlockwise direction, which is taken as the positive direction of the angle.At any time, there is no more than one lug in contact with the terrain when an olig-lugged wheel traverses the reference terrain.Assume θ is equal to α nl /2 at the initial time t 1 when the wheel starts to move.The movement of the olig-lugged wheel rotating around the support lug can be divided into three phases.(1) Wheel center rising phase: this phase starts at moment t 1 , when the support lug starts to be in contact with the terrain, until moment t 2 when the support lug is perpendicular to the terrain.The wheel center moves from the lowest point to the highest point along the y-axis.(2) Wheel center dropping phase: this phase starts at moment t 2 until moment t 3 when the wheel rim starts to be in contact with the terrain.During this phase, the wheel center moves from the highest point to the lowest point along the y-axis.(3) Wheel center height unchanging phase: this phase is from moment t 3 to moment t 4 when the next adjacent lug starts to be in contact with the terrain.The height of the wheel center maintains the lowest invariance.

T2
some lugs contacting terrain surface (SLCT) T3 lugs entering into soil (LEIS) The origin of coordinate system Σo-xy is located at the wheel center when the wheel is stationary.The x-axis is parallel to the horizontal direction, and the y-axis is opposite the gravity direction.Σo-xy does not move with the wheel.The angle velocity of the wheel is ωw, l is the lug height, nl is the number of lugs, and αnl (αnl = 2π/nl) is the angle between two adjacent lugs.
According to the number of lugs on the wheel, the wheel can be divided into two classes: olig-lugged and multi-lugged wheels.An olig-lugged wheel is defined as the lugged wheel whose rim can contact the terrain surface when traversing hard terrain (HT).A multi-lugged wheel, which is often used on rovers, is defined as the wheel equipped with a large number of lugs that keep its rim from making contact with the terrain surface when traversing HT.

Vibration Model of Wheels Moving on Flat Hard Terrain
Flat hard terrain is regarded as the reference terrain.The vertical vibration of a wheel traversing the reference terrain is defined as the vertical reference vibration.
Figure 3 shows the movement phases for a wheel rotating around a supporting lug on a reference terrain.Here, θ (−αnl/2 < θ ≤ αnl/2) is the angle from the direction of the support lug to the gravity direction at the initial time when the wheel starts to move in the anticlockwise direction, which is taken as the positive direction of the angle.At any time, there is no more than one lug in contact with the terrain when an olig-lugged wheel traverses the reference terrain.Assume θ is equal to αnl/2 at the initial time t1 when the wheel starts to move.The movement of the olig-lugged wheel rotating around the support lug can be divided into three phases.(1) Wheel center rising phase: this phase starts at moment t1, when the support lug starts to be in contact with the terrain, until moment t2 when the support lug is perpendicular to the terrain.The wheel center moves from the lowest point to the highest point along the y-axis.(2) Wheel center dropping phase: this phase starts at moment t2 until moment t3 when the wheel rim starts to be in contact with the terrain.During this phase, the wheel center moves from the highest point to the lowest point along the y-axis.(3) Wheel center height unchanging phase: this phase is from moment t3 to moment t4 when the next adjacent lug starts to be in contact with the terrain.The height of the wheel center maintains the lowest invariance.The wheel repeats the above movement phases with the lugs contacting the terrain in turn.The distance h rt (t) between the lowest point of the wheel rim and the reference terrain surface can be expressed as follows: h rt (t) = max(0, (r + l) cos(ϕ rt (t)) − r). ( Sensors 2023, 23, 9752 7 of 32 where r is the wheel radius, l is the lug height, and ϕ rt (t) is the angle between the support lug and gravity direction at the moment t.It can be expressed as The wheel base excitation is q rt (t) = h rt (t) − min(h rt (t)) = (t).Based on Equation (5), the vertical output displacement of the wheel is y rt (t) = q rt (t) = max(0, (r + l) cos(ϕ rt (t)) − r).
For a multi-lugged wheel, the movement of the wheel rotating around the supporting lug can be divided into wheel center rising phase and dropping phase.The wheel center height unchanging phase, in which the wheel rim can be in contact with the terrain, does not exist.The distance h rt (t) is as follows: The output displacement of the wheel is: When the multi-lugged wheel traverses the reference terrain, there are no more than two lugs in contact with the terrain, and the wheel rim cannot be in contact with the terrain surface.Thus, (r + l) cos(ϕ rt (t)) − r is bigger than zero.The h rt (t) calculated using Equation ( 8) is equal to h rt (t) calculated using Equation (5).Equation ( 8) can be regarded as a special case of Equation (5).Whether the lugs on the wheel are olig or multi, the output displacement for a lugged wheel traversing the reference terrain can be expressed uniformly as follows: y rt (t) = q rt (t) = h rt (t) − min(h rt (t)) = max(0, (r + l) cos(ϕ rt (t)) − r) − max(0, cos(α nl /2)(r + l) − r). (10)

Vibration Model of Wheels with All Lugs Contacting Terrain Surface
In an actual situation, the surface of HT is rugged, as shown in Figure 4.The class of wheel-HT interaction is ALCT.The wheel movement on the terrain is similar to the wheel movement on a reference terrain.The characteristics of the wheel movement are irrelative to the state of the wheel at the initial time.Thus, this paper assumes that the lug starts to contact the terrain at the initial moment t 1 .The terrain coordinate system O-XY is located at the terrain surface.The X-axis is parallel to the horizontal direction, and the Y-axis coincides with the y-axis.The roughness function of the terrain is expressed as g(X).The distance h(t) between the lowest point of the wheel rim and the terrain surface is where S x (t) is the wheel displacement along the x-axis, n is the number of times the lugs contact the terrain from t 1 to the moment t, X n is the position of the supporting lug on the x-axis.
where t sn is the moment when the supporting lug starts to depart from the terrain, and t en is the moment when the next lug starts to contact the terrain.
where Sx(t) is the wheel displacement along the x-axis, n is the number of times the lugs contact the terrain from t1 to the moment t, Xn is the position of the supporting lug on the x-axis.
( ) ( ) ( ) where s n t is the moment when the supporting lug starts to depart from the terrain, and e n t is the moment when the next lug starts to contact the terrain.
During the movement of the wheel, the position of the lug adjacent to the supporting lug on the X-axis is ) The distance between the lug adjacent to the supporting lug and the terrain surface is where ( ) is the sum of the distance between the lug adjacent to the supporting lug and wheel center and terrain roughness of position ( ) nl X t on the x-axis.
( ) During the movement of the wheel, the position of the lug adjacent to the supporting lug on the X-axis is The distance between the lug adjacent to the supporting lug and the terrain surface is where h (t) is the sum of the distance between the lug adjacent to the supporting lug and wheel center and terrain roughness of position X nl (t) on the x-axis.Then, n can be expressed as follows: where t + is the moment next to t, and n(0) is equal to 1.The vertical displacement of the wheel is For a multi-lugged wheel, Sensors 2023, 23, 9752 9 of 32 When a multi-lugged wheel traverses HT, since the wheel rim does not contact the terrain, t sn is equal to t en and (r + l) cos(ϕ(t)) + g(X n ) is bigger than g(S x (t)).Thus, the h(t) calculated using Equation ( 24) is equal to h(t) calculated using Equation (11).Equation ( 24) can be regarded as a case of Equation ( 11).Whether the lugs on the wheel are olig or multi, the output displacement for a lugged wheel traversing HT can be expressed uniformly using Equation (22) with h(t) in Equation (11).

Vibration Model of Wheels with Some Lugs Contacting Terrain Surface
In WTIC 2, SLCT, for example, when the wheel moves on gravel terrain (GT), as shown in Figure 5, the lugs will contact gravel during the movement of the wheel.This leads to some of the lugs hanging in the air without contact with the terrain surface or gravel.The angle between two adjacent supporting lugs is no greater than 90 • .The distance between the top point of the i-th lug (lu i ), which will be in contact with the terrain, and the terrain surface is determined as follows: where X (n+i)l (t) is the position of the i-th lug (lu i ), which will be in contact with the terrain on the X-axis.
where t + is the moment next to t, and n(0) is equal to 1.The vertical displacement of the wheel is For a multi-lugged wheel, When a multi-lugged wheel traverses HT, since the wheel rim does not contact the terrain, s n t is equal to e n t and ( ) ( ) x g S t .Thus, the h(t) calculated using Equation ( 24) is equal to h(t) calculated using Equation (11).Equation ( 24) can be regarded as a case of Equation ( 11).Whether the lugs on the wheel are olig or multi, the output displacement for a lugged wheel traversing HT can be expressed uniformly using Equation ( 22) with h(t) in Equation (11).

Vibration Model of Wheels with Some Lugs Contacting Terrain Surface
In WTIC 2, SLCT, for example, when the wheel moves on gravel terrain (GT), as shown in Figure 5, the lugs will contact gravel during the movement of the wheel.This leads to some of the lugs hanging in the air without contact with the terrain surface or gravel.The angle between two adjacent supporting lugs is no greater than 90°.The distance between the top point of the i-th lug (lui), which will be in contact with the terrain, and the terrain surface is determined as follows:  Moreover, the wheel rim may contact the gravel and no lug will be in contact with the terrain some of the time.The minimum distance between the wheel rim and terrain surface is as follows: where S x (t) is the displacement of the wheel center along the X-axis, x w is the position of a point on the wheel rim relative to the wheel center along the x-axis.X (n+i)l (t) and S x (t) can be calculated using the following equations.
. ( 31) x li (t) = (r + l) sin(ϕ(t) + iα nl ) − r sin(ϕ(t) + α w (t)modα nl ). (34) where x wmax0 is the value of x w that maximizes h w (t) at moment t.The number of wheel lugs in contact with the terrain can be expressed as where t + is the moment next to t, h n+i (t) is the sum of the distance between the i-th lug (lu i ) and wheel center and terrain roughness of position X (n+i)l (t) on the x-axis, and m is the number of lugs.The distance h(t + ) is The distance h(0) between the wheel rim and terrain surface at the initial moment t = 0 is equal to max(g(S x (0)), (r + l) cos(ϕ(0)) + g(X 0 )).The vertical displacement of the vibration output is When a multi-lugged wheel moves on GT, some lugs will not be in contact with the surface.Thus, the vertical displacement of the multi-lugged wheel can also be expressed as Equation (41).

Vibration Model of Wheels with Lugs Entering the Soil
In WTIC 3: LEIS, the lugs enter into the soil during wheel movement, for example, when a wheel traverses through sandy terrain (ST), as shown in Figure 6.To facilitate an analysis of the vertical output displacement, assume that a lug starts to enter the soil at the initial moment of the wheel movement.The distance h(t) between the lowest point of the wheel rim and the terrain surface can be expressed as follows: where y l (t) is the position of the top point on the lug starting to enter the soil on the y-axis, β is the angle between the lug starting to enter the soil and gravity direction.The positions of the top points on two adjacent lugs on the y-axis are expressed as y l1 (t) and y l2 (t), respectively.They have the relationship y l2 (t) = y l1 (t − α nl /ω w ) when the wheel-terrain interaction condition, including the terrain properties, slip ratio of the wheel, angular velocity of the wheel, load of the wheel, and so on, remains unchanged.Thus, the position of the top points on the frontmost lug in soil on the y-axis at moment t is y f l (t) = y l (t − INT((ω w t)/α nl ) × α nl /ω w ).It can be used to calculate the distance h(t) at moment t.Therefore, 43) When a multi-lugged wheel moves on GT, some lugs will not be in contact with the surface.Thus, the vertical displacement of the multi-lugged wheel can also be expressed as Equation (41).

Vibration Model of Wheels with Lugs Entering the Soil
In WTIC 3: LEIS, the lugs enter into the soil during wheel movement, for example, when a wheel traverses through sandy terrain (ST), as shown in Figure 6.To facilitate an analysis of the vertical output displacement, assume that a lug starts to enter the soil at the initial moment of the wheel movement.The distance h(t) between the lowest point of the wheel rim and the terrain surface can be expressed as follows: where ( )  According to Equation (43), h(t) is a periodic function with the period α nl /ω w .The value of h(t) (shown in Figure 7a) for an olig-lugged wheel moving on ST varies as follows during the entire process.When the lug starts to contact the terrain at moment t 1 , it decreases until moment t 2 when the lug is parallel to the y-axis and reduces to the minimum.Subsequently, h(t) increases until moment t 3 when the lug departs from the soil and reaches the maximum.Then, it remains constant until the next lug contacts the terrain at moment t 4 .For a multi-lugged wheel, there is more than one lug excavating the soil at the same time.The value of h(t) varies as shown in Figure 7b during the wheel movement.When the lug starts to contact the terrain at moment t 1 , h(t) decreases until moment t 2 when there is one lug in the soil along the y-axis direction and reduces to minimum.Subsequently, h(t) increases until moment t 3 when there are two adjacent lugs symmetrically about the gravity direction and reaches the maximum.Then, it decreases again until the next lug contacts the terrain at moment t 4 .when the wheelterrain interaction condition, including the terrain properties, slip ratio of the wheel, angular velocity of the wheel, load of the wheel, and so on, remains unchanged.Thus, the position of the top points on the frontmost lug in soil on the y-axis at moment t is . It can be used to calculate the distance h(t) at mo- According to Equation (43), h(t) is a periodic function with the period αnl/ωw.The value of h(t) (shown in Figure 7a) for an olig-lugged wheel moving on ST varies as follows during the entire process.When the lug starts to contact the terrain at moment t1, it decreases until moment t2 when the lug is parallel to the y-axis and reduces to the minimum.Subsequently, h(t) increases until moment t3 when the lug departs from the soil and reaches the maximum.Then, it remains constant until the next lug contacts the terrain at moment t4.For a multi-lugged wheel, there is more than one lug excavating the soil at the same time.The value of h(t) varies as shown in Figure 7b during the wheel movement.When the lug starts to contact the terrain at moment t1, h(t) decreases until moment t2 when there is one lug in the soil along the y-axis direction and reduces to minimum.Subsequently, h(t) increases until moment 3 t′ when there are two adjacent lugs symmetrically about the gravity direction and reaches the maximum.Then, it decreases again until the next lug contacts the terrain at moment t4.
(a) (b) The output displacement of the wheel is The output displacement of the wheel is (44)

Simplification of Vibration Outputs for Wheels Traversing Different Terrains
The vibration models built in Equations ( 22), (41), and ( 44) are too complex to analyze the vibration characteristics, thus they need to be simplified.For ALCT, it is known from Equation (20) that the amounts of time each lug contacts the terrain surface is not equal and is influenced by the waves of the terrain.The amplitude of y(t) is slightly larger than that for the wheel traversing the reference terrain surface.However, when the wheel rotates 360 • on T1 and reference terrain with the same angular velocity, the amount of time and the number of times of the lugs contact the terrain are the same.Therefore, the time of each lug contacting the terrain can be regarded as approximately the same.In addition, y(t) can be regarded as a periodic function whose amplitude is variable.Compared with the output y rt (t) of the reference vibration, y(t) can be expressed as follows: where indicates the increase in the amplitude of y(t) caused by the terrain roughness compared to the amplitude of y rt (t).
For SLCT, the amplitude of y(t) is much bigger than y rt (t) for the wheel traversing reference terrain.The number of contact times between wheel and terrain will be less than the number of lugs when the wheel rotates 360 • .So, the frequency of y(t) is lower than n l ω w /2π.The output displacement y(t) of the wheel can be expressed as follows: where C f (t) (C f (t) > 1) indicates the frequency reduction of y(t) incurred by the terrain roughness when compared to the frequency of y rt (t) for the reference terrain.In addition, indicates the increase in the amplitude of y(t) caused by the terrain roughness compared to the amplitude of y rt (t).
For LEIS, every lug can contact the terrain.According to Equations ( 43) and (44), the frequency of y(t) is n l ω w /2π, equal to the frequency of y rt (t).The amplitude of y(t) is smaller than y rt (t) for the wheel traversing reference terrain.Compared with the output y rt (t) of the reference vibration, y(t) can be expressed as where indicates the decrease in the amplitude of y(t) incurred from the lug entering the soil when compared to the amplitude of y rt (t).
Equations ( 45), (46), and (48) show that the amplitude of y(t) is independent of the angular velocity of the wheel.The value of y(t), which can be approximated as a periodic function, can be expressed as follows: where f a is the frequency of output displacement for the vertical vibration, ϕ i is the phase angle of the i-th harmonic, and C y is the constant of y(t) and is equal to −min(h(t)).

Analysis and Extraction of Statistical Vibration Features Independent of the Wheel Velocity
It is difficult to measure vertical displacement of a wheel during movement, whereas vertical acceleration of the wheel is easy to obtain.Therefore, the vibration features of the wheels are extracted from the vertical acceleration signal.The vibration acceleration signal is often described by its statistical features such as mean square value, maximum value, kurtosis, and so on.These statistical features are proportional to wheel speed.If a rover moves on low-unevenness terrain at a high speed, the values of statistical features may be close to that for a rover moving on high-unevenness terrain at a low speed.This can cause low-unevenness terrain to be classified as high-unevenness terrain.In order to solve this problem, new speed-independent features are proposed in this paper.

Extraction of Vibration Features Independent of Wheel Velocity
According to the wheel's vertical vibration displacement y(t) in Equation ( 49), its vertical vibration acceleration can be expressed as a y (t) = ..
The vertical vibration acceleration of the wheel can be regard as an approximate period signal, for example, the vertical acceleration of the wheel is an approximate period signal as shown in Figure 8 when the wheel moves on ST.The Fourier-series expansion of a periodic acceleration signal is described as follows: where C i is the amplitude of the i-th harmonic.

Extraction of Vibration Features Independent of Wheel Velocity
According to the wheel's vertical vibration displacement y(t) in Equation ( 49), its vertical vibration acceleration can be expressed as The vertical vibration acceleration of the wheel can be regard as an approximate period signal, for example, the vertical acceleration of the wheel is an approximate period signal as shown in Figure 8 when the wheel moves on ST.The Fourier-series expansion of a periodic acceleration signal is described as follows: ( ) ( ) where Ci is the amplitude of the i-th harmonic.
According to the analysis in Section II, the changing of wheel-terrain interaction class changes the frequency and amplitude of the vertical vibration acceleration signals of a wheel.The vertical reference acceleration is defined as the vertical acceleration of a wheel when it traverses the reference terrain.Therefore, the vertical reference acceleration can be expressed as Four statistical features are dimensionlessly processed to be independent of wheel speed.They are as follows: (1) The frequency ratio of the vertical vibration acceleration for a wheel traversing an actual terrain and the reference terrain is described as where fa is the frequency of the vertical vibration acceleration for wheels traversing actual terrains with angular velocity ωw, and fart is the frequency of the vertical vibration acceleration for wheels traversing reference terrain with the same angular velocity ωw.
The value of κf is equal to the average value of Cf(t).It is independent of the angular velocity, velocity, and slip ratio of the wheel since Cf(t) is independent of those motion parameters.The value of κf is directly proportional to the frequency of the wheel in contact with the terrain.According to Section II, the frequency of the vibration signal is appropriately equal to the quotient of the number of times the wheel contacts the terrain divided by the cost time when the wheel rotates 360°.According to Equations (45), (46), and (48), κf is equal to 1 for LEIS and ALCT, and κf is smaller than 1 for SLCT.According to the analysis in Section 2, the changing of wheel-terrain interaction class changes the frequency and amplitude of the vertical vibration acceleration signals of a wheel.The vertical reference acceleration is defined as the vertical acceleration of a wheel when it traverses the reference terrain.Therefore, the vertical reference acceleration can be expressed as a yrt (t) = ..
Four statistical features are dimensionlessly processed to be independent of wheel speed.They are as follows: (1) The frequency ratio of the vertical vibration acceleration for a wheel traversing an actual terrain and the reference terrain is described as where f a is the frequency of the vertical vibration acceleration for wheels traversing actual terrains with angular velocity ω w , and f art is the frequency of the vertical vibration acceleration for wheels traversing reference terrain with the same angular velocity ω w .The value of κ f is equal to the average value of C f (t).It is independent of the angular velocity, velocity, and slip ratio of the wheel since C f (t) is independent of those motion parameters.The value of κ f is directly proportional to the frequency of the wheel in contact with the terrain.According to Section 2, the frequency of the vibration signal is appropriately equal to the quotient of the number of times the wheel contacts the terrain divided by the cost time when the wheel rotates 360 • .According to Equations (45), (46), and (48), κ f is equal to 1 for LEIS and ALCT, and κ f is smaller than 1 for SLCT.
(2) The ratio of the sum of all harmonic amplitudes in Equation (50) for a wheel traversing an actual terrain and the reference terrain is described as where A a is the sum of all harmonics' amplitudes of the vertical vibration acceleration for wheels traversing actual terrains with angular velocity ω w , and A art is the sum of all harmonics' amplitudes of the vertical vibration acceleration for wheels traversing reference terrain with the same angular velocity ω w .
The κ A is independent of the angular velocity, velocity, and slip ratio of the wheel because C yi and C yrti are independent of those motion parameters.The κ A is relative to the stiffness and roughness of the terrain.It is directly proportional to the terrain roughness for terrain that does not cause the wheel to sink and inversely proportional to wheel sinkage for terrain that causes the wheel to sink.For ALCT, κ A is slightly greater than 1.For SLCT, κ A is much greater than 1.For LEIS, κ A is less than 1.
(3) The mean value ratio of vertical vibration acceleration for a wheel traversing an actual terrain and the reference terrain is described as follows: where v a is the mean value of the vertical vibration acceleration for wheels traversing actual terrains with angular velocity ω w , and v art is the mean value of the vertical vibration acceleration for wheels traversing reference terrain with the same angular velocity ω w .Due to C y0 and C yrt0 being independent of the wheel's angular velocity, velocity, and slip ratio, the κ mv is independent of those motion parameters.It indicates intensity of the vertical speed change.
(4) The mean square ratio of the vertical vibration acceleration for a wheel traversing an actual terrain and the reference terrain is described as follows: where v a is the mean square value of the vertical vibration acceleration for wheels traversing actual terrains with angular velocity ω w , and v art is the mean square value of the vertical vibration acceleration for wheels traversing reference terrain with the same angular velocity ω w .The value of κ ms is independent of the wheel's angular velocity, velocity, and slip ratio because the numerical value of a y 2 (t)/ω w 4 is independent of those motion parameters.The κ ms is relative to the stiffness and roughness of the terrain, similar to κ A .It can reflect the variation in speed of the signal from the minimum to the maximum value.It is the powerfulness of vibration energy.For ALCT, κ ms is slightly greater than 1.For SLCT, κ ms is much greater than 1.For LEIS, κ ms is smaller than 1.

Analysis of Vibration Features Based on a Single-Wheel Experiment
Experiments were conducted using the wheel-terrain interaction test platform with an accelerometer used to collect the vertical acceleration of a wheel, as shown in Figure 9a.The three motors on the platform are the driving, turning, and dragging motors, and the sensors on the platform include the torque, displacement, and six-axis force/torque sensor.The driving motor drives the wheel forward, and the dragging motor is used to simulate the body movement of the rover.A wheel moving with different slip ratios can be controlled by regulating the two motors at different rotational speeds.The displacement sensor is used to measure wheel sinkage.A counterweight is used to adjust vertical load of the wheel, whereas the actual load is measured using the six-axis force/torque sensor.The radius of the wheel is 140 mm and its width is 150 mm.A wheel with 28 lugs of 15 mm in height is shown in Figure 9b.The vertical load of the wheel is 100 N.For each terrain class, the velocities of the wheel are 10, 14, and 20 mm/s, and the slip ratio of the wheel is 0-0.4.The Mars terrain mainly includes hard, sandy, and gravel surfaces [3].Rovers need to safely traverse those terrains.Motivated by the challenges of rover mobility, HT, GT, ST, hard terrain with gravel (HGT), and sandy terrain with gravel (SGT) were chosen for the experiment, as indicated in Figure 10.HT is made up of hard slate, GT is constructed with a pile of gravel, ST is sand, HGT is formed by placing gravel randomly on slate, and SGT is sand with a small amount of gravel with a random distribution.The wheel-HT interaction pertains to ALCT, wheel-GT/HGT interaction pertains to SLCT, and wheel-ST/SGT interaction belongs to LEIS. Figure 11 shows the extraction results for the vibration features.The values of the vibration features are expressed as a fuzzy set {Small (S), Medium (M), Big (B), Very-Big (VB)}.Table 2 shows a comparison of all vibration features for three WTICs.The vibration features extracted to constitute feature vector P for WTIC recognition are {p1, p2, p3, p4} (p1 = lg(κf), p2 = lg(κA), p3 = lg(κmv), and p4 = lg(κms)) According to P, each WTIC can be distinguished from the others.The Mars terrain mainly includes hard, sandy, and gravel surfaces [3].Rovers need to safely traverse those terrains.Motivated by the challenges of rover mobility, HT, GT, ST, hard terrain with gravel (HGT), and sandy terrain with gravel (SGT) were chosen for the experiment, as indicated in Figure 10.HT is made up of hard slate, GT is constructed with a pile of gravel, ST is sand, HGT is formed by placing gravel randomly on slate, and SGT is sand with a small amount of gravel with a random distribution.The Mars terrain mainly includes hard, sandy, and gravel surfaces [3].Rovers need to safely traverse those terrains.Motivated by the challenges of rover mobility, HT, GT, ST, hard terrain with gravel (HGT), and sandy terrain with gravel (SGT) were chosen for the experiment, as indicated in Figure 10.HT is made up of hard slate, GT is constructed with a pile of gravel, ST is sand, HGT is formed by placing gravel randomly on slate, and SGT is sand with a small amount of gravel with a random distribution.The wheel-HT interaction pertains to ALCT, wheel-GT/HGT interaction pertains to SLCT, and wheel-ST/SGT interaction belongs to LEIS. Figure 11 shows the extraction results for the vibration features.The values of the vibration features are expressed as a fuzzy set {Small (S), Medium (M), Big (B), Very-Big (VB)}.Table 2 shows a comparison of all vibration features for three WTICs.The vibration features extracted to constitute feature vector P for WTIC recognition are {p1, p2, p3, p4} (p1 = lg(κf), p2 = lg(κA), p3 = lg(κmv), and p4 = lg(κms)) According to P, each WTIC can be distinguished from the others.The wheel-HT interaction pertains to ALCT, wheel-GT/HGT interaction pertains to SLCT, and wheel-ST/SGT interaction belongs to LEIS. Figure 11     According to the analysis of Section IIIB, when the WTIC is SLCT, κf is smaller than 1, and κA and κms are greater than 1.The feature vector P for SLCT (such as wheel-HGT/GT interaction) is located in the region {p1, p2, p3, p4 | p1 < lgεfl, lgεAu < p2, lgεmvu < p3, lgεmsu < p4}.When WTIC is LEIS, κf is approximately equal to 1, and κA/κms are smaller than 1.The According to the analysis of Section 3.2, when the WTIC is SLCT, κ f is smaller than 1, and κ A and κ ms are greater than 1.The feature vector P for SLCT (such as wheel-HGT/GT interaction) is located in the region {p 1 , p 2 , p 3 , p 4 | p 1 < lgε fl , lgε Au < p 2 , lgε mvu < p 3 , lgε msu < p 4 }.When WTIC is LEIS, κ f is approximately equal to 1, and κ A /κ ms are smaller than 1.The feature vector P for LEIS (such as wheel-ST/SGT interaction) is in the region {p 1 , p 2 , p 3 , p 4 | lgε fl < p 1 < lgε fu , p 2 < lgε Al , lgε mvl < p 3 <lgε mvu , p 4 < lgε msl }.If P is not located in the above three regions, the corresponding WTIC is classified as an unknown wheel-terrain interaction class (UWTIC).

Analysis of Feature Vector P for the All Lugs Contacting Terrain Surface within the Vicinity of Zero
Owing to the wheels used on the rover all being multi-lugged, the analysis in this section is conducted without considering olig-lugged wheels.When the WTIC is ALCT, the supporting lug of the wheel must depart from the terrain before the adjacent lug along the y-axis.When lug lu 1 contacts the terrain after lug lu 0 along the y-axis and departs from the terrain before the adjacent lug lu 2 along the y-axis, as shown in Figure 13, the vibration displacement caused by lu 1 has the highest amplitude in this case.When lu 1 is along the y-axis, the distance h(t) between the lowest point of the wheel and the terrain is l.The distance h(t) calculated using lu 1 is (r + l) cos(α nl ) − r at the moment when lug lu 2 is along the y-axis.The amplitude of the vibration output displacement caused by rotating around lu 1 is (r + l)(1 − cos(α nl )).Thus, to make every lug have contact with the terrain, the amplitude of the vibration output displacement caused by rotating around each lug is no bigger than (r + l)(1 − cos(α nl )).The amplitude of the reference vibration displacement is (r + l)(1 − cos(α nl /2)).Therefore, C 1 M (t) in Equation ( 45) is no bigger than (1 − cos(α nl ))/(1 − cos(α nl /2)).The following in equalities can be obtained based on Equations ( 52) and (54).( ) Thus, the following equations can be obtained: The parameters εfl, εmvl, εfu, and εmvu are influenced by the noise and accelerometer properties including the sampling frequency and precision.In this study, εfl = 1 − εf, εfu = 1 Thus, the following equations can be obtained: The parameters ε fl , ε mvl , ε fu , and ε mvu are influenced by the noise and accelerometer properties including the sampling frequency and precision.In this study, ε fl = 1 − ε f , ε fu = 1 + ε f , ε f is set as 0.15, ε mvu is set as 10, and the lower limits of ε Al , ε mvl , and ε msl are 0.5.Therefore, the range of the feature vector P for the ALCT is as follows: (59)

Recognition of Wheel-Terrain Interaction Class Based on Vibration Features
According to the value of feature vector P, the WTIC can be recognized.If P satisfies Equation (59), the WTIC is ALCT.According to the analysis in Section 4.1, the value ranges of feature vector P for SLCT and LEIS are expressed as Equation (60) and Equation (61), respectively.
If P is not located within the above three regions, the corresponding WTIC is UWTIC.The programming for WTIC recognition is shown in Figure 14.Table 3 shows the WTIC recognition results for the single-wheel test experiments.The accuracy of the WTIC categorization into the three classes is no less than 96.67%.These results show the ability of the algorithm to distinguish WTICs.When attempting to identify ALCT, the algorithm misclassifies 2 samples of 60 samples as UWTIC, and the recognition accuracy reaches 96.67%.When recognizing SLCT, 2 samples of 120 samples are misclassified as UWTIC, and the recognition accuracy is about 98.34%.Similarly, when recognizing LEIS, 4 samples of 120 samples are misclassified as UWTIC, and the recognition accuracy is approximately 96.67%.The three WTICs are not mutually misclassified.This clearly demonstrates the ability of the algorithm to identify WTICs that induce clearly distinct vibrations on the wheel.
If P is not located within the above three regions, the corresponding WTIC is UWTIC.The programming for WTIC recognition is shown in Figure 14.Table 3 shows the WTIC recognition results for the single-wheel test experiments.The accuracy of the WTIC categorization into the three classes is no less than 96.67%.These results show the ability of the algorithm to distinguish WTICs.When attempting to identify ALCT, the algorithm misclassifies 2 samples of 60 samples as UWTIC, and the recognition accuracy reaches 96.67%.When recognizing SLCT, 2 samples of 120 samples are misclassified as UWTIC, and the recognition accuracy is about 98.34%.Similarly, when recognizing LEIS, 4 samples of 120 samples are misclassified as UWTIC, and the recognition accuracy is approximately 96.67%.The three WTICs are not mutually misclassified.

Collect wheel vibration acceleration
Output WTIC is UWTIC Judge whether P satisfies the Eq. ( 59

Estimation of Terrain Properties Based on the Recognition of Wheel-Terrain Interaction Class
The WTICs can be recognized using the algorithm shown in Figure 14.The terramechanics model is switched corresponding to the recognized WTIC and then used to identify the terrain parameters in this section.

Switchable Terramechanics Model
Yang et al. [19] inferred that the terramechanics models for wheels moving on sand and hard terrain can be similarly expressed, and the parameters of sand and hard terrain can be unified.The terramechanics model equations for wheel-ST interaction can be switched to the equations that have the same form as the Hunt-Crossley contact model by setting the shearing deformation modulus K and soil cohesion c as zero.It can be used for the interaction between a rigid wheel and hard terrain.
The normal force F N , drawbar pull F DP , and resistance moment M R between the wheel and soil can be expressed as follows: where r is the wheel radius; b is the wheel width; z is the wheel sinkage, i.e., deformation of the terrain; c is the soil cohesion; ϕ is the internal friction angle of the soil; N is the sinkage exponent (N = n 0 + n 1 s); k s is expressed as k s = k c /b + k ϕ (where k c is the cohesive modulus of the soil and k ϕ is frictional modulus of soil); and A, B, and C are as follows: where θ 1 (θ 1 = arccos(1-z/r)) is the entrance angle, θ m (θ m = (c 1 + c 2 s)θ 1 ) is the angle of the maximum stress, θ 2 (θ 2 = c 3 θ 1 ) is the departure angle, and c 1 , c 2 , and c 3 are the coefficients of wheel-soil interaction angle.The parameter E is expressed as follows: where K represents the shearing deformation modulus, and j is the shearing displacement.When the wheel contacts hard terrain, the terrain has no shearing deformation.Thus, K is equal to zero.In addition, c can be set to zero according to Yang et al. [19].Equations ( 62)-( 64) are as follows: For interaction between the wheel and hard terrain, θ 1 tends toward zero, and thus the coefficient of friction µ tends toward tan ϕ.The parameter λ c is the coefficient used to balance the driving torque.When the parameters K and c are appropriately assigned, the terramechanics model can be used for a rover traversing hard terrain.Therefore, the terramechanics model for wheel-ST interaction can be switched to be suitable for wheel-HT interaction.

Unified Identification Model for Terrain Parameters
Let X = tan ϕE, Y = rbk s z N , and Q = rbcE; the following equation can be obtained: where W and T D are the wheel load and driving torque of the motor, respectively.Equations ( 73)-( 75) are a decoupled terramechanics model regarding unknown terrain parameters K s , n 0 , n 1 , c 1 , c 2 , c, ϕ, and K, which can be used for terrain parameter identification.

Identification of Terrain Parameters Based on Terramechanics Model Switching
Table 4 shows the correspondence between WTICs and setting values of terrain parameters for terramechanics model switching.WTIC recognition is used to judge whether some parameters should be set to zero to achieve terramechanics model switching.If the WTIC is ALCT, c is set to zero, K is set to zero, E is set to 1.This is because there is no shear deformation or backward accumulation of the soil for ALCT.GT is a typical terrain pertaining to T2 (wheel-T2 interaction is SLCT).When the wheel traverses GT, gravel will move slightly backward and accumulate, which is far less than the accumulation degree of sand.Thus, if the WTIC is SLCT, according to [19], c is set to zero.The sinkage and slip ratio of the wheel are necessary parameters for identification of the terrain parameters in the terramechanics model.Wheel sinkage can be detected using machine vision methods or through a kinematic analysis of the rover suspension [21,22].When WTIC is not LEIS, sinkage is set as 0.01mm for terrain parameter identification.The wheel slip ratio can be estimated using a machine vision method [23].The vertical load W, driving torque T D , and drawbar pull F DP can be measured by installing a six-axis force/torque sensor on the wheel [18].
According to Equation (73), F DP is a function of W, T D , θ 1 , s, and the unknown parameters c 1 and c 2 .Here, c 1 and c 2 can be identified according to the measured F DP , W, T D , θ 1 , and s.According to Equation (74), W is a function of T D , θ 1 , z, s, identified parameters c 1 and c 2 , and unknown parameters K s , n 0 , and n 1 .The three parameters K s , n 0 , and n 1 can be identified using the measured W, T D , θ 1 , z, s, and identified c 1 , c 2 .According to Equation (75), T D is a function of W, θ 1 , s, identified parameters c 1 and c 2 , and unknown parameters c, ϕ, and K. Parameters c, ϕ, and K can be identified using the measured W, T D , θ 1 , s, and identified values of c 1 and c 2 .
The least squares method is used to achieve identification for unknown terrain parameters.When identifying c 1 and c 2 , xe denotes the parameter vector {c 1 , c 1 } to be identified; xdata denotes the input data vector {s, F N , T D }; ydata denotes the measured F DP ; m denotes the length of xdata and ydata; and F denotes the function of Equation (73).When identifying K s , n 0 , and n 1 , xe denotes the parameter vector {K s , n 0 , n 1 } to be identified; xdata denotes the input data vector {s, T D }; ydata denotes the measured F N ; and F denotes the function of Equation ( 74).When identifying c, ϕ, and K, xe denotes the parameter vector {c, ϕ, K} to be identified; xdata denotes the input data vector {s, F N }; ydata denotes the measured T D ; and F denotes the function of Equation (75).The ultimate goal is to find the vector xe that best fits Equation (76): Figure 15 shows the programming of the terrain parameter identification for different types of terrain.

Experiments for Wheel-Terrain Interaction Class Recognition and Terrain Parameter Identification Using a Planetary Rover Prototype
The method proposed in this paper is experimentally validated using a six-wheeled Planetary Rover Prototype.The collected vibration acceleration of the wheel is used to recognize the WTIC.WTIC recognition guides the rover to switch the terramechanics model for identifying terrain parameters in the model using the algorithm shown in Figure 15.The results for the WTIC recognition and identification of terrain parameters in the terramechanics model are discussed in this section.

Introduction of Planetary Rover Prototype and Terrain Environment
The Planetary Rover Prototype is shown in Figure 16.An accelerator mounted on the front left leg of the rover, near the joint with the wheel axle, is used to collect the vertical acceleration of the wheel to study the vibration-based WTIC recognition.A camera is used to monitor the robot's wheel and detect the wheel sinkage and slip ratio.The camera is

Experiments for Wheel-Terrain Interaction Class Recognition and Terrain Parameter Identification Using a Planetary Rover Prototype
The method proposed in this paper is experimentally validated using a six-wheeled Planetary Rover Prototype.The collected vibration acceleration of the wheel is used to recognize the WTIC.WTIC recognition guides the rover to switch the terramechanics model for identifying terrain parameters in the model using the algorithm shown in Figure 15.The results for the WTIC recognition and identification of terrain parameters in the terramechanics model are discussed in this section.

Introduction of Planetary Rover Prototype and Terrain Environment
The Planetary Rover Prototype is shown in Figure 16.An accelerator mounted on the front left leg of the rover, near the joint with the wheel axle, is used to collect the vertical  Three different surfaces were selected for the experiments: GT, HT (tiled bricks), and ST, as shown in Figure 17.HT pertains to terrain T1, GT pertain T2, and ST pertains to terrain T3.As such, they represent a wide range of hard soft to hard, and a roughness ranging from flat to rugged, thus better illustratin all accuracy of the methods presented.

Recognition of Wheel-Terrain Interaction Class for Planetary Rover Prototype
On each of the terrain surfaces, the rover drove back and forth thrice a velocity of 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, and 0.01 m/s.Durin movement, 32-bit vibration acceleration samples were collected using the shown in Figure 14 at a frequency of 128 kHz.The signal is processed by dow and low-pass filtering to extract vibration features.Figures 18 and 19 show th signal samples for three classes of terrain.Once the single-terrain datasets wer multi-terrain datasets were collected with the rover driving from ST to HT to gle traversal, both forward and reverse.
A single dataset consists of vibration data from an about 55 s traversal Three different surfaces were selected for the experiments: GT, HT (tiled with stone bricks), and ST, as shown in Figure 17.HT pertains to terrain T1, GT pertains to terrain T2, and ST pertains to terrain T3.As such, they represent a wide range of hardnesses from soft to hard, and a roughness ranging from flat to rugged, thus better illustrating the overall accuracy of the methods presented.
ensors 2023, 23, x FOR PEER REVIEW rotational DOFs, installed on the robot body.A gyroscope is attached to detect the rotation angle of the camera.The position of the camera can be the relation of the arm and robot body.Every wheel is equipped w force/torque sensor to detect the drawbar pull force FDP, vertical load W (eq force FN), lateral force FL, driving moment TD, steering moment MS, and ov ment MO.The radius of the wheel is 155 mm and its width is 150 mm.Th lugs of 10 mm in height.Three different surfaces were selected for the experiments: GT, HT (ti bricks), and ST, as shown in Figure 17.HT pertains to terrain T1, GT pert T2, and ST pertains to terrain T3.As such, they represent a wide range of ha soft to hard, and a roughness ranging from flat to rugged, thus better illustr all accuracy of the methods presented.

Recognition of Wheel-Terrain Interaction Class for Planetary Rover Prototyp
On each of the terrain surfaces, the rover drove back and forth thric velocity of 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, and 0.01 m/s.Du movement, 32-bit vibration acceleration samples were collected using t shown in Figure 14 at a frequency of 128 kHz.The signal is processed by d and low-pass filtering to extract vibration features.Figures 18 and 19 show signal samples for three classes of terrain.Once the single-terrain datasets w multi-terrain datasets were collected with the rover driving from ST to HT gle traversal, both forward and reverse.
A single dataset consists of vibration data from an about 55 s travers acceleration signals are segmented with 15 s long intervals to recognize

Recognition of Wheel-Terrain Interaction Class for Planetary Rover Prototype
On each of the terrain surfaces, the rover drove back and forth thrice at a constant velocity of 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, and 0.01 m/s.During the rover movement, 32-bit vibration acceleration samples were collected using the accelerator shown in Figure 14 at a frequency of 128 kHz.The signal is processed by down-sampling and low-pass filtering to extract vibration features.Figures 18 and 19 show the vibration signal samples for three classes of terrain.Once the single-terrain datasets were collected, multiterrain datasets were collected with the rover driving from ST to HT to GT in a single traversal, both forward and reverse.
A single dataset consists of vibration data from an about 55 s traversal period.The acceleration signals are segmented with 15 s long intervals to recognize WTIC.Table 5 show the results for the WTIC recognition of single-terrain datasets.They demonstrate an excellent performance of the WTIC recognition method proposed in this paper in distinguishing the WTICs.The results show that 133 samples of 144 wheel-HT interaction (which belongs to ALCT) are recognized as ALCT, and the recognition accuracy reaches 92.36%.Ten samples of wheel-GT interaction (which belongs to SLCT) are misclassified as UWTIC, i.e., 93.06% of the wheel-GT interactions are classified as SLCT.All wheel-ST interaction samples are recognized as LEIS, and the accuracy attains 100%.No misidentification occurred among the three WTICs.The results show that the WTIC recognition algorithm has a high level of accuracy and can effectively distinguish the WTICs.Figure 20 shows information regarding the vibration acceleration and WTIC recognition results for a distinct single traversal of all three terrain classes shown in Figure 15.The rover is traveling from an area of ST to an area of HT and then on to an area of GT.The data are also divided into 15 s interval segments for WTIC recognition.During the entire movement, the Planetary Rover Prototype operates for a total of 315 s, thus, the acceleration is divided into 21 segments for WTIC recognition.The segmentation is shown in the dotted square in Figure 18, and the WTIC recognition results are represented by the mark at the bottom of the block.These results clearly show the effectiveness of the algorithm in recognizing multiple WTICs during a single traversal.This figure also illustrates the most common error mechanism of the algorithm, which occurs when the rover transitions between homogenous terrain regions.Figure 20 shows information regarding the vibration acceleration and WTIC recognition results for a distinct single traversal of all three terrain classes shown in Figure 15.The rover is traveling from an area of ST to an area of HT and then on to an area of GT.The data are also divided into 15 s interval segments for WTIC recognition.During the entire movement, the Planetary Rover Prototype operates for a total of 315 s, thus, the acceleration is divided into 21 segments for WTIC recognition.The segmentation is shown in the dotted square in Figure 18, and the WTIC recognition results are represented by the mark at the bottom of the block.These results clearly show the effectiveness of the algorithm in recognizing multiple WTICs during a single traversal.This figure also illustrates the most common error mechanism of the algorithm, which occurs when the rover transitions between homogenous terrain regions.Figure 20 shows information regarding the vibration acceleration and WTIC recognition results for a distinct single traversal of all three terrain classes shown in Figure 15.The rover is traveling from an area of ST to an area of HT and then on to an area of GT.The data are also divided into 15 s interval segments for WTIC recognition.During the entire movement, the Planetary Rover Prototype operates for a total of 315 s, thus, the acceleration is divided into 21 segments for WTIC recognition.The segmentation is shown in the dotted square in Figure 18, and the WTIC recognition results are represented by the mark at the bottom of the block.These results clearly show the effectiveness of the algorithm in recognizing multiple WTICs during a single traversal.This figure also illustrates the most common error mechanism of the algorithm, which occurs when the rover transitions between homogenous terrain regions.

Terrain Parameter Identification for Planetary Rover Prototype
The wheel sinkage and slip ratio are necessary parameters for terrain parameter identification and can be detected using machine vision methods shown in our earlier studies [22,23].

Introduction of Sinkage Detection for Wheels Moving on Soft Terrain
Ref. [22] detailed the detection method of wheel sinkage, as shown in Figure 21.It proposes a new definition of wheel sinkage for rough terrain.Wheel sinkage has been defined as a vector using the wheel-soil boundary.A new vision-based method for sinkage detection is proposed in the paper.The saturation of wheel-terrain image is adjusted using a dynamic piecewise non-linear adjustment to enhance the color contrast between the wheel and soil regions.Through a binarizing process, the wheel section in the color-enhanced image is segmented to acquire an image in which the wheel region is white and other regions are black.Thus, the outline of the wheel region can be extracted from the binary image using the edge detection method.The wheel section outline is corrected to represent the wheel rim operating as a circle.The distance between points of the wheel rim and wheel center is equal to the wheel radius.And the distance between points of the wheel-soil boundary and wheel center is less than the wheel radius.The wheel-soil boundary is extracted based on this characteristic.The sinkage is calculated using the wheel-soil boundary based on the sinkage definition.

Terrain Parameter Identification for Planetary Rover Prototype
The wheel sinkage and slip ratio are necessary parameters for terrain parameter identification and can be detected using machine vision methods shown in our earlier studies [22,23].

Introduction of Sinkage Detection for Wheels Moving on Soft Terrain
Ref. [22] detailed the detection method of wheel sinkage, as shown in Figure 21.

Terrain Parameter Identification for Planetary Rover Prototype
The wheel sinkage and slip ratio are necessary parameters for terrain parameter identification and can be detected using machine vision methods shown in our earlier studies [22,23].

Introduction of Sinkage Detection for Wheels Moving on Soft Terrain
Ref. [22] detailed the detection method of wheel sinkage, as shown in Figure 21.It proposes a new definition of wheel sinkage for rough terrain.Wheel sinkage has been defined as a vector using the wheel-soil boundary.A new vision-based method for sinkage detection is proposed in the paper.The saturation of wheel-terrain image is adjusted using a dynamic piecewise non-linear adjustment to enhance the color contrast between the wheel and soil regions.Through a binarizing process, the wheel section in the color-enhanced image is segmented to acquire an image in which the wheel region is white and other regions are black.Thus, the outline of the wheel region can be extracted from the binary image using the edge detection method.The wheel section outline is corrected to represent the wheel rim operating as a circle.The distance between points of the wheel rim and wheel center is equal to the wheel radius.And the distance between points of the wheel-soil boundary and wheel center is less than the wheel radius.The wheel-soil boundary is extracted based on this characteristic.The sinkage is calculated using the wheel-soil boundary based on the sinkage definition.It proposes a new definition of wheel sinkage for rough terrain.Wheel sinkage has been defined as a vector using the wheel-soil boundary.A new vision-based method for sinkage detection is proposed in the paper.The saturation of wheel-terrain image is adjusted using a dynamic piecewise non-linear adjustment to enhance the color contrast between the wheel and soil regions.Through a binarizing process, the wheel section in the color-enhanced image is segmented to acquire an image in which the wheel region is white and other regions are black.Thus, the outline of the wheel region can be extracted from the binary image using the edge detection method.The wheel section outline is corrected to represent the wheel rim operating as a circle.The distance between points of the wheel rim and wheel center is equal to the wheel radius.And the distance between points of the wheel-soil boundary and wheel center is less than the wheel radius.The wheel-soil boundary is extracted based on this characteristic.The sinkage is calculated using the wheel-soil boundary based on the sinkage definition.

Introduction of Wheel Slip Ratio Detection
Ref. [23] detailed a detection method for wheel slip ratio shown in Figure 22.The wheel slip ratio was estimated using the same wheel-terrain images as used in the sinkage detection.At different moments, the wheel position in the wheel-terrain image is the same, while the position of the wheel-trace boundary is different.Thus, the wheel speed can be estimated from the displacement of the wheel trace relative to the wheel at different moments.Based on this idea, Lv et al. [23] built an estimation model for the slip ratio, which is expressed as: where ∆S is the displacement of the wheel, ∆u is the average displacement of the mark points, and r 1 is the radius of the points.Ref. [23] detailed a detection method for wheel slip ratio shown in Figure 22.The wheel slip ratio was estimated using the same wheel-terrain images as used in the sinkage detection.At different moments, the wheel position in the wheel-terrain image is the same, while the position of the wheel-trace boundary is different.Thus, the wheel speed can be estimated from the displacement of the wheel trace relative to the wheel at different moments.Based on this idea, Lv et al. [23] built an estimation model for the slip ratio, which is expressed as: ( ) 1 arccos 1 2 where ΔS is the displacement of the wheel, Δu is the average displacement of the mark points, and r1 is the radius of the points.As shown in Figure 22a, the boundaries of the wheel trace were extracted from wheel-terrain images acquired at adjacent moments and were matched to estimate the forward displacement ΔS of the wheel.Meanwhile, the mark points (in this paper, the screws on the wheel are regarded as mark points) of the wheel were also extracted from wheel-terrain images acquired at adjacent moments and were matched to estimate the average displacement Δu of the mark points relative to the wheel center.The slip ratio can then be calculated using Equation (77) according to ΔS and Δu.
When the rover is traversing on HT or GT, there is no obvious wheel trace on the terrain, but there are SURFF feature points that do not exist on ST.The wheel trace boundary extraction and matching proposed by Lv et al. [23] for slip ratio estimation are replaced by the terrain SURFF feature point extraction and matching.The forward displacement ΔS of the wheel can be estimated by detecting the SURFF feature points in the wheelterrain images acquired at two adjacent moments, thereby estimating the wheel slip ratio, as shown in Figure 22b.As shown in Figure 22a, the boundaries of the wheel trace were extracted from wheelterrain images acquired at adjacent moments and were matched to estimate the forward displacement ∆S of the wheel.Meanwhile, the mark points (in this paper, the screws on the wheel are regarded as mark points) of the wheel were also extracted from wheelterrain images acquired at adjacent moments and were matched to estimate the average displacement ∆u of the mark points relative to the wheel center.The slip ratio can then be calculated using Equation (77) according to ∆S and ∆u.
When the rover is traversing on HT or GT, there is no obvious wheel trace on the terrain, but there are SURFF feature points that do not exist on ST.The wheel trace boundary extraction and matching proposed by Lv et al. [23] for slip ratio estimation are replaced by the terrain SURFF feature point extraction and matching.The forward displacement ∆S of the wheel can be estimated by detecting the SURFF feature points in the wheel-terrain images acquired at two adjacent moments, thereby estimating the wheel slip ratio, as shown in Figure 22b.

Identification Results of Terrain Parameters for Three Classes of Terrain
The slip ratio of the object wheel can be adjusted by varying the speed of the object wheel and maintaining other wheels at a constant speed.For each terrain class, the ratio is controlled at 0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, and 0.6.The actual slip ratio is detected visually.
Terrain parameters can be divided into three groups, namely, P I = {c 1 , c 2 }, P II = {k s , n 0 , n 1 }, and P III = {c, ϕ, K}, which are directly related to angles of θ m and θ 2 , normal stress, and shearing stress, respectively.The results of terrain parameter identification using the algorithm in Figure 15 are shown in Table 6.Figures 23 and 24 show comparison of the parameters P II = {k s , n 0 , n 1 } and P III = {c, ϕ, K} for three terrain classes, respectively, to clearly illustrate the difference in terrain parameters between the three terrain classes.The slip ratio of the object wheel can be adjusted by varying the speed of the object wheel and maintaining other wheels at a constant speed.For each terrain class, the ratio is controlled at 0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, and 0.6.The actual slip ratio is detected visually.
Terrain parameters can be divided into three groups, namely, PI = {c1, c2}, PII = {ks, n0, n1}, and PIII = {c, φ, K}, which are directly related to angles of θm and θ2, normal stress, and shearing stress, respectively.The results of terrain parameter identification using the algorithm in Figure 15 are shown in Table 6.Figures 23 and 24 show comparison of the parameters PII = {ks, n0, n1} and PIII = {c, φ, K} for three terrain classes, respectively, to clearly illustrate the difference in terrain parameters between the three terrain classes.Wheel-terrain interaction forces are predicted for the wheels without considering the rover dynamics or force coordination among the wheels during the rover movement to simplify the calculation of them for easy control.The estimated drawbar pull force FDP, wheel load W, and driving torque TD with identified terrain parameters are compared with the sensor-measured results as shown in Figures 25-27.Their estimated relative errors are shown in Table 7. From Figures 25-27 and Table 7, the predicted wheel-terrain interaction forces are close to the measured values.Thus, the simplified predicted wheelterrain force can be used for optimal control (such as to maximize the wheel traction or minimize the power consumption) for rovers.The slip ratio of the object wheel can be adjusted by varying the speed of the object wheel and maintaining other wheels at a constant speed.For each terrain class, the ratio is controlled at 0, 0.05, 0.1, 0.15, 0.2, 0.3, 0.4, 0.5, and 0.6.The actual slip ratio is detected visually.
Terrain parameters can be divided into three groups, namely, PI = {c1, c2}, PII = {ks, n0, n1}, and PIII = {c, φ, K}, which are directly related to angles of θm and θ2, normal stress, and shearing stress, respectively.The results of terrain parameter identification using the algorithm in Figure 15 are shown in Table 6.Figures 23 and 24 show comparison of the parameters PII = {ks, n0, n1} and PIII = {c, φ, K} for three terrain classes, respectively, to clearly illustrate the difference in terrain parameters between the three terrain classes.Wheel-terrain interaction forces are predicted for the wheels without considering the rover dynamics or force coordination among the wheels during the rover movement to simplify the calculation of them for easy control.The estimated drawbar pull force FDP, wheel load W, and driving torque TD with identified terrain parameters are compared with the sensor-measured results as shown in Figures 25-27.Their estimated relative errors are shown in Table 7. From Figures 25-27 and Table 7, the predicted wheel-terrain interaction forces are close to the measured values.Thus, the simplified predicted wheelterrain force can be used for optimal control (such as to maximize the wheel traction or minimize the power consumption) for rovers.Wheel-terrain interaction forces are predicted for the wheels without considering the rover dynamics or force coordination among the wheels during the rover movement to simplify the calculation of them for easy control.The estimated drawbar pull force F DP , wheel load W, and driving torque T D with identified terrain parameters are compared with the sensor-measured results as shown in Figures 25-27.Their estimated relative errors are shown in Table 7. From Figures 25-27 and Table 7, the predicted wheel-terrain interaction forces are close to the measured values.Thus, the simplified predicted wheel-terrain force can be used for optimal control (such as to maximize the wheel traction or minimize the power consumption) for rovers.6 shows that the k s of HT is much larger than that of ST because the internal bonding force of HT is much greater than that of ST and has a large rigidity; thus, the terrain does not undergo substantial deformation when a wheel moves on it.There is no sinkage of the wheel.The sinkage exponents n 0 and n 1 have different characteristics for three classes of terrain.The n 0 of HT is close to 0.5, whereas n 1 is approximately zero.Thus, the sinkage exponent N is not affected by wheel slip ratio.The normal force model of the wheel-HT interaction is a non-linear elastic model as shown in Equation (67).The tangential bond of HT is extremely large.During the interaction between the wheel and HT, the soil will not be damaged.There is no shear deformation or backward accumulation for the terrain.Therefore, there is no shear stress in the force of the wheel-terrain interaction, and only normal stress occurs.From Figure 25a, the fitting curve of estimated wheel load W using the identified parameters is extremely close to the fitting curve of the measured data.The relative error of estimated wheel load W is no more than 11.68% as shown in Table 7.
The drawbar pull force can be calculated using the friction model as shown in Equation (68).When the wheel slips, the drawbar pull force of the wheel is equal to the friction between the wheel and the terrain.The parameters c and K are both zero, ϕ represents the external friction angle of the terrain, and tanϕ is the friction coefficient between the terrain and wheel.
Because the wheel sinkage is close to zero, the maximum stress angle θ m is calculated as equal to zero according to the coefficients of the wheel-soil interaction angle.
Figure 25a shows that the estimated values of F DP and T D differ greatly from the measured values when slip ratio is 0. And the maximum relative error of the estimated drawbar pull force is 65.8%, as shown in Table 7.However, the estimated F DP and T D are close to measured data except when the slip ratio is 0, as shown in Figure 25b.The relative errors of estimated drawbar pull force and driving torque in Figure 25b are both no more than 15.52%.When the slip ratio is zero, there is no relative movement between the wheel and terrain, and the drawbar pull force (friction force) of the wheel is static friction force.At this time, the actual drawbar pull force (friction force) is smaller than the dynamic friction force calculated according to ϕ (the friction coefficient between wheel and terrain is tan(ϕ)).Thus, the errors of the estimated drawbar pull force and the driving torque are large when slip ratio is zero.Figure 25b shows a comparison among the estimated values of the drawbar pull force, the normal load, and driving torque and the sensor measurement values when the slip ratio is bigger than zero.The fitting curve of estimated data is close to the fitting curve of measured data.The relative errors of estimated wheel-HT interaction force with identified terrain parameters are less than 16% when slip ratio is not zero.

Analysis of Parameter Identification Results of Gravel Terrain
Table 6 shows that GT's k s is close to HT's, GT has a large rigidity, and the terrain does not undergo substantial deformation when a wheel moves on it.GT's n 0 is close to 0.5 and n 1 is about 0.054.The value of n 1 is relatively small, but the value of k s is very large.The normal force is sensitive to the sinkage exponent.Therefore, the normal force in Figure 26 has obvious variation.
Under the action of wheels, gravel will slightly move backward and accumulate, which is far less than the accumulation degree of sand.There is shear displacement between pieces of gravel and shear stress in the force of the terrain on the wheel.Therefore, K is not 0 in the shear parameter, and ϕ represents the external friction angle of the terrain.
Since the wheel sinkage is close to 0, the maximum stress angle calculated from coefficients of the wheel-soil interaction angle is approximately equal to 0.
The drawbar pull force, normal load, and driving torque of the wheel are estimated using Equations ( 62)-(64) with the identified parameters.Figure 26 shows that the fitting curves of estimated wheel interaction force are close to the fitting curves of measured data.The maximum errors of estimated drawbar pull force, wheel load, and driving torque are 11.85%, 6.67%, and 8.46% respectively.

Analysis for Parameter Identification Results of Sandy Terrain
Table 6 shows that k s of ST is much smaller than that of GT/HT.The bonding force between sand grains is small, and thus the stiffness of sand is low.The wheel easily sinks into the soil when the rover traverses over ST.The n 0 of ST is about 0.8169 and n 1 is about 1.2957.Thus, the sinkage exponent N is greatly affected by the wheel slip ratio.The normal stress model is a variable exponential non-linear elastic model as shown in Equation (62).
During the interaction between the wheel and terrain, the soil will be displaced, causing the compression under the wheel and accumulation behind the wheel.The terrain places a large amount of shear stress on the wheel, which is not symmetrically distributed along the y-direction of the wheel.The shear stress is also an important factor in calculating the supporting force of the terrain applied to the wheel.
Owing to the obvious shear deformation, both c and K are greater than zero and much larger than those of GT.In addition, ϕ is the internal friction angle of the soil.
From Figure 27, it can be seen that the relative error of the estimated drawbar pull force, wheel load, and driving torque using Equations ( 62)-( 64) with the identified parameters is small, and the fitting curves of the estimated data are extremely close to the fitting curves of the measured data.Table 7 shows that the maximum relative errors of the estimated drawbar pull force, wheel load, and driving torque are 8.11%, 7.64%, and 7.75%, respectively.
According to the above analysis of experiment results, it can be seen that the relative errors of estimated wheel-terrain interaction force with identified terrain parameters are less than 16%, 12%, and 9% for rovers traversing hard, gravel, and sandy terrain, respectively.

Discussion and Conclusions
This paper has been focused on speed-independent vibration-based WTIC recognition.The terramechanics model was switched for terrain parameter identification according to the recognized WTIC.Concluding remarks of this study are provided as follows: (1) Three vibration models of rigid lugged wheels have been built and analyzed for three WTICs.The three vertical vibration displacements are all approximately periodic functions.Therefore, the vibration acceleration of a rigid lugged wheel was approximated as a periodic signal.Its frequency, amplitude, mean value, and mean square value were extracted as vibration features describing the wheel-terrain interaction response characteristics and processed to make them non-dimensional and independent of the wheel speed.
(2) A vibration-based method for WTIC recognition without requirement of training data is proposed.The regions of the vibration feature space for three WTICs are analyzed.The WTICs are determined by judging which feature space region the wheel-terrain interaction is located in.Experiment results of the single-wheel test show the simple and effective method has recognized 292 samples out of 300 samples correctly, so the recognition accuracy is about 97%.
(3) The proposed WTIC recognition method can help rovers to switch the terramechanics model when traversing hard or soft terrain by modulating the parameters c and K appropriately.Terrain parameters in the switched terramechanics model are identified using the least squares method, and model parameters have been successfully identified for various terrain classes, particularly for hard terrain/rock and soft sandy terrains that are similar to the main terrain types of Mars.The experiment results obtained using a Planetary Rover Prototype moving on the three terrain classes have demonstrated that the WTIC recognition, terramechanics model switching, and parameter identification are effective.The relative errors of estimated wheel-terrain interaction force with identified terrain parameters are less than 16%, 12%, and 9% for rovers traversing hard, gravel, and sandy terrain, respectively.
A vibration-based method for WTIC recognition is proposed in this paper.Previous planetary terrain analysis studies usually identify terrain types, such as ST, GT, HT, SGT, HGT, and so on [4][5][6][7][8].When rovers move on some of these terrain types (such as ST and SGT), the WTIC is the same.It is difficult to give comprehensive and accurate statistics on the correspondence between terrain type and WTIC.Thus, these methods cannot recognize WTIC, select terramechanics models, and identify terrain parameters.In addition, terrain classification may result in misidentification between different types of terrain.It can lead to incorrect selection of terramechanics models.The WTIC recognition method proposed in this paper does not involve misidentification between different WTICs.It is more conducive to the selection of terramechanics models and the identification of terrain parameters.
In terms of terrain parameter identification, a method for identifying the terrain parameters when the rover moves on various terrains is proposed in this paper.Previous studies on terrain parameter identification (such as [16][17][18]) focused on identifying the terrain parameters of ST without identifying the terrain parameters of other types of terrain.They have not been able to switch and select terramechanics models.If the terramechanics model of wheel-ST interaction is used for terrain parameter identification of other types of terrain, it can lead to incorrect terrain parameter identification results.In terms of identifying the terrain parameters of ST, the method proposed in this paper is similar to [18].The accuracy of identified terrain parameters is basically the same, and the relative errors of estimated wheel-terrain interaction force with identified terrain parameters are also close.Compared to previous studies on terrain parameter identification, the proposed method can achieve the terrain parameter identification of various types of terrain.It has improved the parameter identification ability.
As future work, the identification of terrain parameters in a terramechanics model in real time for planetary rovers traversing various types of terrain will be further studied.The parameters identified using this method will be used to improve the estimation accuracy of wheel-terrain interaction force in the control model for a planetary rover.And they can be used to calculate the wheel-driving torque for evaluating the terrain traversability [24].In addition, to improve the autonomous mobility of a rover, a terrain-environment perception system of a planetary rover can be built by combining this study with vision-based terrain classification.The system can help a planetary rover construct a knowledge base of the visual terrain features and terramechanics properties and roughly estimate the terrain parameters' value ranges in the terramechanics model for terrains that the rover will move on.

Figure 1 .
Figure 1.Model of a rigid lugged wheel traversing on any terrain.
frequency of the undamped free mass-spring-damper system.According to [20], λ is less than 0.75.

Figure 1 .
Figure 1.Model of a rigid lugged wheel traversing on any terrain.

Figure 3 .
Figure 3. Movement phases of a wheel rotating around the supporting lug when the wheel traverses flat hard terrain.

Figure 4 .
Figure 4. Wheel motion states of rotating around the supporting lug on natural HT.

Figure 4 .
Figure 4. Wheel motion states of rotating around the supporting lug on natural HT.

Figure 5 .
Figure 5. Wheel motion states of rotating around the supporting lug on GT.Figure 5. Wheel motion states of rotating around the supporting lug on GT.

Figure 5 .
Figure 5. Wheel motion states of rotating around the supporting lug on GT.Figure 5. Wheel motion states of rotating around the supporting lug on GT.

Figure 6 .
Figure 6.Wheel motion states when moving on sandy terrain.

(
y t is the position of the top point on the lug starting to enter the soil on the yaxis, β is the angle between the lug starting to enter the soil and gravity direction.The positions of the top points on two adjacent lugs on the y-axis are expressed as 1 ( ) l y t and

Figure 6 .
Figure 6.Wheel motion states when moving on sandy terrain.

Figure 7 .
Figure 7.The distance h(t) between wheel rim's lowest point and terrain surface, (a) when an oliglugged wheel moves on ST, (b) when a multi-lugged wheel moves on ST.The black solid line represents the acceleration signal, the red dashed line represents the time point, and the black dashed line represents the coordinate grid.

Figure 7 .
Figure 7.The distance h(t) between wheel rim's lowest point and terrain surface, (a) when an olig-lugged wheel moves on ST, (b) when a multi-lugged wheel moves on ST.The black solid line represents the acceleration signal, the red dashed line represents the time point, and the black dashed line represents the coordinate grid.

Figure 8 .
Figure 8. Sample of vibration acceleration when a wheel moves on ST.

2 )Figure 8 .
Figure 8. Sample of vibration acceleration when a wheel moves on ST.

Sensors 2023 ,Figure 9 .
Figure 9. Experimental equipment, (a) is the wheel-soil interaction test platform and (b) is the experimental wheel.

Figure 9 .
Figure 9. Experimental equipment, (a) is the wheel-soil interaction test platform and (b) is the experimental wheel.

Sensors 2023 ,Figure 9 .
Figure 9. Experimental equipment, (a) is the wheel-soil interaction test platform and (b) is the experimental wheel.

35 Figure 11 .
Figure 11.Values of vibration features for different terrain classes.Vibration features are expressed as the fuzzy set {Small (S), Medium (M), Big (B), Very-Big (VB)}.

4 Figure 11 .
Figure 11.Values of vibration features for different terrain classes.Vibration features are expressed as the fuzzy set {Small (S), Medium (M), Big (B), Very-Big (VB)}.

Figure 12 .
Figure 12.The regions of feature space formed by the vibration features for different WTICs, with (a) the projection in plane p1p2, (b) the projection in plane p1p3, (c) the projection in plane p1p4, (d) the projection in plane p2p3, (e) the projection in plane p2p4, (f) the projection in plane p3p4.

Figure 12 .
Figure 12.The regions of feature space formed by the vibration features for different WTICs, with (a) the projection in plane p 1 p 2 , (b) the projection in plane p 1 p 3 , (c) the projection in plane p 1 p 4 , (d) the projection in plane p 2 p 3 , (e) the projection in plane p 2 p 4 , (f) the projection in plane p 3 p 4 .

κ A ≤ ( 1 −Figure 13 .
Figure 13.The wheel-terrain interaction case in which vibration displacement caused by lu1 has the biggest amplitude.

Figure 13 .
Figure 13.The wheel-terrain interaction case in which vibration displacement caused by lu1 has the biggest amplitude.

Figure 14 .
Figure 14.Flow chart for recognition of wheel-terrain interaction class.

Figure 15 .
Figure 15.Flow chart for identifying parameters of various types of terrain.

Figure 15 .
Figure 15.Flow chart for identifying parameters of various types of terrain.
wheel to study the vibration-based WTIC recognition.A camera is used to monitor the robot's wheel and detect the wheel sinkage and slip ratio.The camera is mounted on a supporting arm, with three moving degrees of freedom (DOFs) and three rotational DOFs, installed on the robot body.A gyroscope is attached to the camera to detect the rotation angle of the camera.The position of the camera can be acquired with the relation of the arm and robot body.Every wheel is equipped with a six-axis force/torque sensor to detect the drawbar pull force F DP , vertical load W (equal to normal force F N ), lateral force F L , driving moment T D , steering moment M S , and overturning moment M O .The radius of the wheel is 155 mm and its width is 150 mm.The wheel has 24 lugs of 10 mm in height.Sensors 2023, 23, x FOR PEER REVIEW rotational DOFs, installed on the robot body.A gyroscope is attached to the detect the rotation angle of the camera.The position of the camera can be acq the relation of the arm and robot body.Every wheel is equipped with force/torque sensor to detect the drawbar pull force FDP, vertical load W (equa force FN), lateral force FL, driving moment TD, steering moment MS, and overt ment MO.The radius of the wheel is 155 mm and its width is 150 mm.The w lugs of 10 mm in height.

Sensors 2023 ,Figure 20 .
Figure 20.Vibration acceleration and terrain classification results for distinct single traversals of all three terrain classes.

Figure 20 .
Figure 20.Vibration acceleration and terrain classification results for distinct single traversals of all three terrain classes.

Sensors 2023 ,Figure 20 .
Figure 20.Vibration acceleration and terrain classification results for distinct single traversals of all three terrain classes.

Figure 22 .
Figure 22.Flow charts of wheel slip ratio detection, with (a) for the wheel traversing ST, and (b) for the wheel traversing GT/HT.

Figure 22 .
Figure 22.Flow charts of wheel slip ratio detection, with (a) for the wheel traversing ST, and (b) for the wheel traversing GT/HT.

Figure 24 .
Figure 24.Comparison of parameters c, φ, and K for three terrains.

Figure 23 .
Figure 23.Comparison of parameters k s , n 0 , and n 1 for three terrains.

Figure 24 .
Figure 24.Comparison of parameters c, φ, and K for three terrains.

Figure 24 .
Figure 24.Comparison of parameters c, ϕ, and K for three terrains.

Figure 26 .
Figure 26.Wheel-terrain interaction force for the Planetary Rover Prototype moving on GT.

Figure 27 .Figure 25 .Figure 25 .
Figure 27.Wheel-terrain interaction force for the Planetary Rover Prototype moving on ST.

Figure 26 .
Figure 26.Wheel-terrain interaction force for the Planetary Rover Prototype moving on GT.

Figure 27 .Figure 26 .Figure 25 .
Figure 27.Wheel-terrain interaction force for the Planetary Rover Prototype moving on ST.

Figure 26 .
Figure 26.Wheel-terrain interaction force for the Planetary Rover Prototype moving on GT.

Figure 27 .Figure 27 .
Figure 27.Wheel-terrain interaction force for the Planetary Rover Prototype moving on ST.

Table 1 .
Correspondence between terrain classes and the WTICs.

Table 2 .
The comparison of all vibration features for different WTICs.

Table 3 .
WTIC recognition results of 300 samples for the single-wheel test experiments.

Table 3 .
WTIC recognition results of 300 samples for the single-wheel test experiments.

Table 4 .
Correspondence between terrain class and setting values of terrain parameters for terramechanics model switching.

Table 5 .
WTIC recognition results for Planetary Rover Prototype test experiments.

Table 5 .
WTIC recognition results for Planetary Rover Prototype test experiments.

Table 5 .
WTIC recognition results for Planetary Rover Prototype test experiments.

Table 6 .
Results of terrain parameter identification for Planetary Rover Prototype test experiments.

Table 6 .
Results of terrain parameter identification for Planetary Rover Prototype test experiments.

Table 6 .
Results of terrain parameter identification for Planetary Rover Prototype test experiments.

Table 7 .
Relative error of data fitting for F DP , W, and T D .