## Abstract

This article presents a new method of detecting incipient immobilization for a wheeled mobile robot operating in deformable terrain with high spatial variability. This approach uses proprioceptive sensor data from a four-wheeled, rigid chassis rover operating in poorly bonded, compressible snow to develop canonic, dynamical system models of robot’s operation. These serve as hypotheses in a multiple model estimation algorithm used to predict the robot’s mobility in real time. This prediction method eliminates the need for choosing an empirical wheel–terrain interaction model, determining terramechanics parameter values, or for collecting large training datasets needed for machine learning classification. When tested on field data, this new method warns of decreased mobility an average of 1.8 m and 2.9 s before the rover is completely immobilized. This system also proves to be a reliable predictor of immobilization when evaluated in simulated scenarios of rovers with passive suspension maneuvering in more variable terrain.

## 1 Introduction

A primary design consideration in developing any ground vehicle is the terrain it is expected to encounter and traverse. This constraint needs to be balanced with requirements for the power budget, efficiency of operation, cost, and complexity. For rovers meant to operate on deformable terrain, such as sand or snow, these requirements result in an intractable optimization problem due to the wide range of terrain parameters that are likely. It is impossible to ensure that the rover’s mechanical design alone allows it to gain traction in all types of terrain encountered. This renders the rover continually susceptible to immobilization.

For missions where an autonomous vehicle needs to operate in remote regions, an immobilization event could necessitate an expensive, potentially dangerous, and resource-intensive recovery of the robot, or if a rescue is not possible, it could result in the end of the scientific survey being conducted and loss of the vehicle. The application on which this article is focused is the autonomous collection of ground-based measurements on the ice sheets of Greenland and Antarctica. Measurements of snow temperature, albedo, and mass balance in these remote areas are necessary to build climate models and predict sea-level rise [1,2]. While these data can largely be collected via remote sensing, the high accuracy of ground-based measurements are necessary to calibrate and validate the data from satellite or airborne sensors. Utilizing autonomous rovers to traverse the ice sheets and collect these ground-truth data provides higher spatial resolution than in situ point measurements, such as automatic weather stations. They also have the potential to replace manned traverses, which are often limited in scope by high costs or dangerous conditions. These concepts also extend to autonomous operation of ground vehicles on other worlds, where immobilization can abruptly end a mission.

Given the advantages of autonomous surveys in these remote areas, we have designed and deployed three lightweight rovers designed to tow scientific instruments on snow. All three rovers follow the same basic scheme, with four independently driven wheels instead of tracks to minimize cost, weight, and complexity and to maximize efficiency. With each design iteration, the range of terrains over which the rover could sustain mobility is expanded. The first rover, *Cool Robot*, has a rigid chassis and off-the-shelf ATV tires as the tractive elements. The rigid chassis means that the rover can lose traction on undulating terrain [3]. *Yeti* incorporated a central pivot that allows for rotation of the front and rear axles about the roll axis enabling it to maintain four-wheel contact even when traversing sastrugi in Antarctica, but Yeti can still get mired in soft or low cohesion snow [4]. *FrostyBoy* incorporates wider, custom wheels in order to reduce ground pressure, limit sinkage, and improve trafficability in this type of snow [5]. While these rovers have proved largely reliable in executing surveys on the order of tens to hundreds of kilometers in Antarctica and Greenland [6–9], their intermittent immobilizations prevent them from being a viable option for longer duration scientific traverses of the ice sheets. Despite the modification of the rover’s tractive elements, patches of low stiffness or low cohesion terrain can put an end to a traverse if the robot does not adjust to the changing conditions.

This article presents a new method by which proprioceptive sensor data can be used to detect incipient immobilization. With enough warning of an impending loss of mobility, the rover’s control algorithm could adapt to optimize traction, choose an alternate route, or temporarily release a towed load to reduce required drawbar pull, making reliable detection of impending immobilization a need for robots driving on soft terrain. Prior methods of incipient immobilization detection focus on one of the two approaches: (1) estimation of terrain parameters associated with the semi-empirical Wong–Reece model [10] and subsequently assessing trafficability of the terrain given these parameter estimates, and (2) machine learning approaches that, absent any dynamic model of the vehicle or terrain, use signatures of a sensor suite to detect impending immobilization. Both methods have drawbacks, as detailed in Sec. 2. The former lacks mathematical uniqueness and has computational complexity and sensing requirements rendering it difficult for real-time implementation. The latter requires ample training data for both good and poor mobility conditions, which presents challenges in obtaining such data. This article presents a new method by which proprioceptive sensor data can be used to detect incipient immobilization. The method incorporates simple, canonical models identified from a modest data set, as predictors of the near-term evolution of robot motion. By using these models, this new approach incorporates Bayesian decision-making to determine which model is most likely to be governing the dynamics in a given instance. To develop the method, we use field data from immobilization events and simulated data from a truth model derived from these field data. We present data from field tests used to develop the model that show the complex dynamics of immobilization of a four-wheeled vehicle, whereas previous studies have considered immobilization characteristics through controlled, single-wheel studies. Our study shows that models of immobilization require attention to the characteristics of the full four-wheel vehicle dynamics interacting with variable terrain. The truth model is used to assess the algorithm in immobilization scenarios expected as vehicle design is altered, but that were not directly encountered in the field. Section 2 presents the relevant background in terramechanics theory and the pertinent examples of previous efforts to detect impending mobility loss in real time. Section 3 describes how the immobilization data were collected and details the development of the truth model. Section 4 outlines all aspects of the incipient immobilization detection algorithm (IDeA) and discusses the effectiveness of this algorithm when applied to both the collected and simulated data. Finally, Sec. 5 summarizes the impact of this research and the further studies that need to be conducted to validate this method for a range of terrain types and rover designs.

## 2 Background

*W*is the weight of the robot,

*F*

_{x}is the drawbar pull,

*T*

_{r}is the resistive torque, and

*σ*(

*θ*) and

*τ*(

*θ*) are the pressure distribution and shear distribution below the wheel, respectively (see Fig. 1) [10]. The physical dimensions of the wheel are captured with

*r*

_{w}being the wheel radius and

*b*, the wheel width. Equation (1) defines the amount of sinkage the robot wheels will experience by balancing the force exerted on the wheels by the terrain with the normal force of the robot. With this sinkage known, the shear force that can be applied by the wheel to the terrain before failure can be calculated (first term of Eq. (2)). This force must be greater than the compaction resistance force in opposition to the robot’s motion (second term of Eq. (2)), in order for drawbar pull to be positive and for the robot to successfully traverse the terrain. The resistive torque (3) is the minimum torque that must be transmitted to the wheels in order to overcome the terrain forces.

*k*, defines a linear regime where pressure is proportional to sinkage, but a sharp increase in pressure occurs as the snow compacts against a much denser layer within the snowpack at

*z*

_{m}[17]. This general paradigm makes sense when thinking about how a layer of snow, exposed to an extended period of sunlight and wind, could form a sintered, dense crust over which a series of precipitation events deposit fresh, unconsolidated snow. Fully defining terrain also requires a shear ring or shear plate test to determine terrain’s cohesion,

*c*, angle of internal friction,

*ϕ*, and the shear displacement modulus,

*K*

_{s}, to determine the shear stress at failure,

*τ*, according to the Mohr–Coulomb failure criterion:

*j*is the shear displacement [16]. Yet, this still does not result in the full suite of parameters necessary to fully define the wheel–terrain interaction. The pressure distribution below the wheel,

*σ*as a function of

*θ*, remains unknown without direct measurement using a spatially resolved pressure sensor or estimates of the angle of maximum pressure,

*θ*

_{m}, and the angle at which the wheel loses contact with the terrain,

*θ*

_{2}(see Fig. 1). Vehicle-mounted bevameters and vehicle-mounted single-wheel field testers have been developed to evaluate these properties semi-autonomously in situ [18–20]. A similar method that removed the need for additional, specialized equipment on-board the vehicle was developed by NASA in using one of the Mars rover Sojourner’s wheels as a trenching device to determine a subset of the terrain parameters [21]. While these methods prove to be effective in characterizing terrain and comparing actual mobility with predicted mobility from the model, they cannot prevent immobilization. These measurements require the vehicle to stop, losing all forward momentum, and assume low spatial variability of the terrain since they take measurements only at discrete locations. Since snow has been found to have a high degree of spatial variability from bevameter tests performed manually [22], this method would not work in predicting mobility of the rover even one wheel rotation beyond the test location.

*θ*

_{2}to 0 deg, and presumes

*θ*

_{m}to be exactly halfway between

*θ*

_{1}and

*θ*

_{2}. Wong and Reece proved experimentally that this symmetry only occurs at one value of slip, with empirically determined parameters

*c*

_{1}and

*c*

_{2}defining how the angle of maximum stress,

*θ*

_{m}moves based on the wheel’s slip ratio,

*i*, according to the equation [10]:

*i*, defined as follows:

*r*

_{w}being the wheel radius,

*ω*being the wheel speed, and

*v*being the longitudinal velocity of the vehicle. Since all terrain–wheel interaction forces are based on these stress distributions, approximations in terrain models can have significant impacts on determining mobility. This method also is non-unique in determining terrain parameters using sensor measurements and from the need for a significant number of data points before the model converges. Even with a perfectly accurate understanding of the normal pressure distribution below a wheel, as well as an exact calculation of drawbar pull for that wheel, both the cohesion parameter,

*c*, and the angle of internal friction,

*ϕ*, contribute to the shear stress according to Eq. (5), each with significantly different impacts on mobility at varying slip ratios. Therefore, multiple data points need to be employed to extract each parameter’s respective value. In addition, torque and normal load are typically not directly measured for wheeled mobile robots (WMRs), as this would require a multi-axis force sensor in each wheel adding significant complexity, cost, and mass to a rover. Since these inputs are estimated from a quasi-static analysis based on mass distribution and from current supplied to the motor, the approximations will be inherently noisy or inaccurate. In trying to account for sensor noise, Iagnemma et al. corrupted the inputs with white noise, and this, combined with the issue of non-uniqueness, necessitated at least 30 measurements before the model converged [23]. Depending on the sampling rate and the speed of the rover, the need for 30 samples before convergence could require a long stretch of homogeneous terrain to obtain an accurate estimate of mobility and provide an ample time window for the robot to lose all traction. In Ref. [23], the robot was traveling at 5–10 cm/s, the sensor sampling rate was 5 Hz, and the terrain can be assumed to be somewhat homogeneous, resulting in this method providing adequate warning. However, for

*FrostyBoy*traveling at 1–2 m/s, sampling at 10 Hz, and on highly variable terrain, all four wheels could be immobilized by the time the model converges.

Other online estimation methods have relied on specific maneuvers to be performed by the rover before parameters can be estimated. Specifically, Ojeda et al. performed turned at varying yaw rates while measuring motor currents in order to develop an analog to the shear–stress versus shear–displacement curve [24]. Another study conducted by Liu et al. [25] extracts terrain parameters only when a WMR is skidding (negative slip) down a slope. Both these methods prove limited in their application since requiring the rover to perform turns could incite an immobilization event, and downhill slopes are not guaranteed (nor likely on a flat ice sheet), and the parameters estimated during these maneuvers may not be applicable to upcoming terrain.

These methods have the same four drawbacks in predicting mobility:

They require an assumption that the terrain is homogeneous to make accurate predictions about future mobility.

They calculate a subset of the terrain parameters, with average values assumed for the remaining terrain parameters.

They require simplifications of the pressure and shear stress distributions at the wheel–terrain interface to solve the Wong–Reece equilibrium equations in real time.

The terrain parameters are inextricably linked in their contribution to mobility and cannot be uniquely identified from single force measurements or estimates.

To solve these issues and predict mobility in nonhomogeneous terrain, Ray [15] employed Bayesian multiple model estimation (MME) to evaluate the most likely set of terrain parameters that could describe the robot’s current and past sensor readings. In this method, Ray uses multiple hypotheses of almost a complete set of terrain parameters (*θ*_{2} is assumed to be 0 deg) as inputs to the Wong–Reece equilibrium equations to generate expected forces from the terrain on a rigid-wheeled, lightweight rover [15]. The forces estimated through an extended Kalman–Bucy filter from the proprioceptive sensor measurements of the wheels can then be compared to those generated from the hypotheses, and the most likely set of terrain parameters can be determined based on their conditional probability given the measurements. Ray [15] tested this method on simulations of a constant terrain type and found that this method converged within two iterations of the MME algorithm to the closest terrain parameter hypothesis even when no hypothesis was an exact match of the terrain parameters being simulated. As Cook [26] found, this rapid convergence meant Ray’s method was effective in detecting terrain parameters quickly in nonhomogeneous snow, allowing traction control methods to be developed to maintain mobility in changing conditions. Cook added the shear deformation modulus, *K*_{s}, in his terrain parameter hypotheses, meaning one fewer parameter that is fixed to the average value found in the literature [26]. The third drawback of previous methods is also solved, since the force–slip curves generated for each hypothesis can be calculated numerically ahead of deployment, meaning that no simplifications to the Wong–Reece equations are necessary for the algorithm to run in real time. However, the issue of non-uniqueness still plagues this method as Ray [15] found overlap in differentiating high-cohesion soils in her simulations, and Cook’s algorithm [26] was computationally expensive as it ran 320 parallel Kalman filters, one for each of the hypotheses in a large set of terrain parameter combinations.

While the MME method is effective at rapidly converging on a set of terrain parameters, the forces these terrain parameters predict and the resulting mobility judgment are still rooted in the Wong–Reece terramechanics model. This model was developed for heavier vehicles, with larger wheels, operating in sands and soils. Modifications to the model have been proposed to better fit the data from lightweight vehicle tests, but these often add parameters to the existing set and are not specific to snow, and none of these adaptations have been universally adopted. In addition, observable effects of the vehicle–terrain dynamics that cause immobilization are not considered in these models.

The primary effect that is not considered explicitly in the methods presented earlier is the dynamic sinkage that occurs from excessive wheel slip and the resultant excavation of terrain. The model relies on an accurate definition of the pressure distribution below a rigid wheel to accurately predict a vehicle’s mobility. This distribution defines the overall sinkage and therefore the compaction resistance and resistive torque applied to the wheel by the terrain. However, while Reece’s Eq. (6) takes into account the radial movement of the angle of maximum stress, *θ*_{m}, based on the wheel’s slip, there is no equation that redefines the pressure distribution after slip sinkage occurs. The boundaries of the contact patch are also poorly understood when slip sinkage occurs. Even in Ray’s study [15], *θ*_{2} was assumed to be directly below the wheel axle, but it has been experimentally shown in single-wheel studies that at higher slip ratios, the backward flowing zone of soil under the wheel enlarges until at 100% slip, and there is no forward flowing zone and the contact patch extends well behind the point directly below the axle (*θ*_{2} is much less than 0 deg) [16]. This slip sinkage effect was the primary cause of *Cool Robot*’s immobilizations observed in the 2017 field season near Summit Station, Greenland, and is shown in Fig. 2. The robot would encounter a patch of low cohesion, low stiffness snow. While making very little forward progress, the robot wheels would slip and dig themselves deeper until the chassis was touching the snow surface. While simple models predicting additional sinkage due to slip as proportional to initial sinkage have been proposed and serve as a good approximation, it does not consider how slip sinkage develops over time. The vehicle–terrain interaction has been observed to be more complicated when this slip sinkage effect occurs [27].

With so many deviations of terrain models from observation of how lightweight, wheeled vehicles operate in snow, a new immobilization detection method needs to be developed that will capture the dynamics of the system without relying on unknown and marginally applicable terrain parameters.

One such method of predicting mobility avoids the problems inherent in measuring or estimating terrain parameters by instead using a machine learning approach. Trautmann and Ray [28] used classified proprioceptive sensor data during normal operation and immobilization events of a rover on snow to train a machine learning algorithm. Instead of a terramechanics model, the method relies on supervised learning of a support vector machine for deciding whether the terrain being traversed is a threat to mobility. This method was proven experimentally using *Yeti*, instrumented with encoders and current sensors on each motor, an IMU, and an optical ground speed sensor [28]. The machine learning approach proved to be very effective in predicting a loss of traction; however, the large training data set required renders it impractical. This method is “expensive” for the operator, as Trautmann hand labeled the robot sensor measurements to characterize them to be one of the three conditions: fully mobile, almost immobilized, or immobilized. In addition, field seasons to collect these data are typically short and costly, and many examples of the robot’s performance on both immobilizing and traversable terrains are required for a robust model.

There is a need for a mobility prediction method that does not rely on estimating terrain parameters, can run in real time, and is not predicated on an enormous amount of training data. It also must be reliable enough to detect an immobilization in a wide variety of terrain conditions, as snow has been observed to be highly variable. The method proposed here addresses this need by using the limited experimental data set to create low-dimensional models of the dynamics of a vehicle under both incipient immobilization and under normal operation in order to leverage a Bayesian approach to mobility prediction that runs in real time due to low computational complexity. The method avoids the need for a large amount of training data for a machine learning approach as shown in Ref. [28] and also avoids parameter estimation approaches that require complex, nonlinear dynamics models as given in Refs. [15,26].

## 3 Data Collection and Truth Model

The first step in developing an algorithm to detect an incipient immobilization is to collect proprioceptive sensor data from naturally evolving immobilizations using an instrumented rover, similar to the in situ testing performed by Trautmann and Ray [28]. These data are used both to develop a truth model and also to evaluate the immobilization detection algorithm.

The truth model is used to provide simulation data to evaluate the immobilization detection algorithm, supplementing the field data and accommodating robots with geometries that differ from the one robot used to acquire field data. The truth model is rooted in Wong– Reece terramechanics but enhanced by real-world observations to mimic observed dynamics. The truth model developed here models dynamics observed in field data but that are not elicited by Wong–Reece terramechanics alone and that highlight the importance of the vehicle dynamics in concert with terramechanics in inducing immobilization evolution observed in the field experiments.

Both the collected data of the rover losing traction and the simulated data of the modeled rover becoming immobilized are used to develop the IDeA—a detection algorithm that can warn the rover of an impending mobility loss. This proposed method has many advantages when compared to the state-of-the-art immobilization prediction approaches. The models for nominal operation and those characteristic of an immobilization are derived directly from data measured with the pertinent robot operating in the range of terrain types for which it was designed. These canonical models produce an immobilization signature and are free from all underlying assumptions regarding the underlying terramechanics or the relationship between terrain parameters and immobilization. This ensures the immobilization detection algorithm is model agnostic, which is particularly important to lightweight rovers in low cohesion terrain, such as snow, where terramechanics model parameters have wide variations. The method also goes beyond single-wheel studies, by considering the collective role of all four wheels in a loss of traction.

With an accurate truth model, these algorithms can also be tested on new scenarios and new terrains that are not tested with the instrumented robot. This allows for a small data set, meaning a less costly, shorter field season for collecting these data, to provide a broadly applicable incipient immobilization detection method.

### 3.1 Data Collection.

Measurements of robot ground speed, roll, pitch, and yaw rates; motor currents; motor speeds; longitudinal and lateral accelerations; as well as drawbar load were logged during outdoor tests on relevant terrains. The robot *FrostyBoy* [29] was used for this testing. It is similar to *Cool Robot* [3] and *Yeti* [4] in that it has four electric motors that are powered by lithium ion batteries, optionally charged by solar panels towed behind on a sled. *FrostyBoy* is heavier than the other two robots with a mass of about 90 kg. Its chassis is constructed of T-channel aluminum extrusion, with a pivot on the front axle to give the robot some compliance about the roll axis and ±9 deg about the yaw axis to reduce lateral bulldozing resistance in a skid-steer turn. The heavier chassis is necessary to support the wider wheels. This robot has been operationally tested in Greenland during a 2018 and a 2019 field season [5]. After the 2019 field season, the robot was outfitted with sensors to measure the parameters listed earlier. The ground speed, as well as accelerations of the robot, are derived using the VectorNav VN-200 GPS-aided Inertial Navigation System mounted on the robot chassis with built-in Kalman filtering algorithms to improve accuracy. This device can measure accelerations with a resolution less than 0.5 mg and has velocity accuracy of less than 0.05 m/s. Roll and pitch accuracies have a standard deviation of 0.03 deg, and heading accuracy, when moving, has a standard deviation of 0.2 deg. The motor currents and motor speeds were read using the RoboteQ FBL2360 brushless DC motor controller. An INTERFACE STA-3 S-beam load cell with a rated load of 4900 N and an accuracy of ±1 N was mounted on the tow point of *FrostyBoy* and measured the instantaneous force required to move a weighted sled. All sensor data were logged at a sampling rate of 10 Hz by a Campbell Scientific CR1000X configured with a CFM100 CompactFlash Memory Module. Final sinkage of each wheel can be inferred from pictures taken after an immobilization.

The instrumented robot was tested in two locations outside Lebel-sur-Quévillon in Québec, Canada, in loose, unconsolidated snow that measured more than 55 cm in depth. The goal of these tests was to gather sensor measurements during an immobilization event. At each test site, *FrostyBoy* was run in an open-loop configuration in which all four motor controllers sent the same voltage, around 40% of the rated voltage, to each of the four motors. The rover was left to operate at this voltage, tracking relatively straight, until it became naturally immobilized by what can only be deduced are natural variations in the snow properties, since the snow appeared highly uniform and completely undisturbed from the surface. To better characterize the terrain, four pressure-sinkage curves were generated from load versus distance data logged by a custom device designed to drive a flat puck of a known area down into the snow. Two penetration tests were performed at locations where the rover was able to maintain mobility, and the other two were taken next to the front wheels after *FrostyBoy* lost traction.

Important findings from these data were the easily discernible patterns and signatures from slip, wheel speed, and torque measurements that have the potential to predict and prevent immobilization. The rover became immobilized frequently at both test sites, typically traveling less than 20 m before becoming stuck, even without any towed load. The pattern observed in every one of these immobilization events was an increase in torque and decrease in wheel speed for two of the wheels and a decrease in torque and increase in wheel speed for the other two wheels. The wheels showing the same behavior were always at a diagonal to one another. An example is shown in Fig. 3, with the front right and rear left wheels exhibiting the same behavior, while the front left and rear right wheels exhibit the opposite behavior. As shown in this example, two of the wheels’ torques climb to the motor controllers’ current limit restriction as their rotational speeds drop to zero (the oscillations observed in the torque record is an artifact of the motor controllers’ programming to prevent the motor from overheating by dropping the current by a set amount each time the current limit of 40 amps is reached). The other two wheels’ rotational speeds increase as their torques decrease, appearing to almost reach a no-load condition that also results in their slip ratios jumping up quickly to 100%. Figures 4 and 5 show observed characteristics of these immobilization events. One is that the two wheels with high slip ratios are observed to excavate the terrain below the wheels, evidenced by the buildup of snow behind them. It is also evident from these photos that there is significant snow “trapped” inside the frames of the two stalled wheels, which adds to the resistance torque. The robot pitching up is another signature of this event, which can be seen in Fig. 5. While this added resistance, due to gravity, can be accounted for using the pitch angle measured by the IMU and knowledge of the robot’s center of mass, it serves as valuable insight into modeling the dynamics of these immobilizations. Specifically, this example shows that immobilization events are not only impacted by the vehicle–terrain interaction but also depend on vehicle’s rigid-body dynamics. The cause of immobilization can therefore not be completely understood by single-wheel studies, and full system testing must be performed to be able to detect incipient immobilization, which is an important conclusion of this study.

To further understand the cause of this complete loss of mobility, we examine the data collected from the penetration tests to better grasp the pressure distributions below the wheels. The four sets of pressure versus sinkage measurements are shown in Fig. 6 with the best-fit Preston–Thomas (4) pressure-sinkage curves plotted over them. It is evident from these data that the immobilizing snow (Data 3 and 4) is highly compressible, having low values of *k*, with deeper denser layers, and higher values of *z*_{m}, than the terrain on which the rover maintained mobility (Data 1 and 2). Figure 6 also shows the variability of the snow, with values of *k* ranging from 19.2 to 65.0 kPa/m and values of *z*_{m} ranging from 0.29 to 0.54 m, with a difference of ≈10 kPa/m between Data 3 and 4, taken just 2 m apart. The variation in snow conditions further reinforces how robust an incipient immobilization detection method must be if it is to ensure rover mobility across deep snow.

### 3.2 Truth Model

#### 3.2.1 Dynamics.

While the collected data are invaluable in developing an incipient immobilization detection algorithm, it is limited in its scope. Immobilization events were only recorded during 2 days of testing in Mar. 2020, with two types of snow and two sets of environmental conditions. This lack of control over the independent variables influencing the rover’s mobility is the nature of in situ testing. Also, a key goal of this work is to generate a mobility prediction method without requiring the exhaustive training data sets necessary for machine learning. Therefore, a truth model is developed to complement the data collected from this relatively small sample size of conditions. Once this model, rooted in the terramechanics models of Bekker, Wong, and Reece, has been adapted to incorporate the dynamics seen in the real-world data, the prediction method can then be tested using simulations of the rover. The truth model allows for more extensive testing of the reliability of incipient immobilization detection for conditions not encountered during field testing, leading to a more comprehensive incipient immobilization detector.

*FrostyBoy*. These body-fixed equations of motion are given by

*F*

_{x}is drawbar pull force on the four wheels in the longitudinal,

*x*, direction and the

*T*

_{r}is the resistive torques for all four wheels applied by the terrain. The geometry and properties of the robot are captured in the terms

*m*, being the robot’s mass of 90 kg,

*t*

_{w}, being the track width of the vehicle (see Fig. 7),

*I*

_{zz}, the yaw moment of inertia about its center of mass, and

*I*

_{ω}, the wheel moment of inertia about its center of rotation. The lateral forces on the wheels are not captured explicitly, but a restoring moment,

*M*

_{res}, is applied to the robot based on its yaw rate, $\psi \u02d9$, to account for the lateral resistance from the terrain on the sidewalls of the wheels. Also, since the rover’s yaw is controlled by skid steering rather than by changing the heading of the wheels or axles relative to the body of the robot, the lateral velocity,

*v*

_{y}, is not included in the model. The

*B*

_{ω}term takes into account the mechanical damping in the drivetrain. The values for these parameters measured from

*FrostyBoy*and used in the truth model are outlined in Table 1.

FrostyBoy parameter values | |||
---|---|---|---|

Description | Variable | Value | Units |

Mass | m | 90 | kg |

Wheel radius | r_{w} | 0.265 | m |

Track width | t_{w} | 1.142 | m |

Distance from COM to front axle | d_{f} | 0.579 | m |

Distance from COM to rear axle | d_{r} | 0.515 | m |

Height of COM above axles | h | 0.140 | m |

Rover moment of inertia about z-axis at COM | I_{zz} | 36.05 | kg m^{2} |

Wheel moment of inertia about axis of rotation | I_{ω} | 0.415 | kg m^{2} |

Restoring moment constant | M_{res} | 0.3 | Nm/(rad/s) |

Mechanical damping constant | B_{ω} | 0.327 | Nm/(rad/s) |

FrostyBoy parameter values | |||
---|---|---|---|

Description | Variable | Value | Units |

Mass | m | 90 | kg |

Wheel radius | r_{w} | 0.265 | m |

Track width | t_{w} | 1.142 | m |

Distance from COM to front axle | d_{f} | 0.579 | m |

Distance from COM to rear axle | d_{r} | 0.515 | m |

Height of COM above axles | h | 0.140 | m |

Rover moment of inertia about z-axis at COM | I_{zz} | 36.05 | kg m^{2} |

Wheel moment of inertia about axis of rotation | I_{ω} | 0.415 | kg m^{2} |

Restoring moment constant | M_{res} | 0.3 | Nm/(rad/s) |

Mechanical damping constant | B_{ω} | 0.327 | Nm/(rad/s) |

*ρ*, to be 0 since this observed variable in the test data had a mean of zero and an absolute maximum less than 1.5 deg. Also, this maximum roll angle would result in less than 1.5 N of weight shift between the left and right wheels on a given axle considering the geometry and weight of

*FrostyBoy*at 0 deg pitch. Therefore, for a given axle, the minimum sinkage among the two wheels is applied to both wheels in the next time-step. A pitch angle is calculated using these minimum sinkage values for each axle according to the following equations:

*z*values are the calculated sinkages and

*d*

_{f}and

*d*

_{r}are the distances from the center of mass to the front and rear axles, respectively. The distribution of weight,

*F*

_{zf}and

*F*

_{zr}, between the two axles is then calculated by taking pitch and the effect of acceleration into account:

*z*

_{fl}>

*z*

_{fr}in the previous time-step, then the normal force from the terrain,

*F*

_{zfl}, for the current time-step is calculated using the given minimum sinkage,

*z*

_{f}, which is equal to

*z*

_{fr}. Then, the normal force supporting the right, front wheel,

*F*

_{zfr}, is calculated from

*z*

_{fr}, being iteratively calculated according to the balance of forces Eq. (1). These calculations are critical in defining the wheel–terrain interaction that corresponds to observation, as they serve as the basis for the reaction forces from the terrain.

#### 3.2.2 Terrain.

The terrain resistance torques and forces used in the truth model are calculated from the Wong–Reece terramechanics model. While the drawbacks and limitations of these empirical equations in simulating a lightweight vehicle in snow are discussed in Sec. 2, it serves as a widely accepted, standard foundation for an investigation of mobility. Therefore, the range of values for the empirical parameters that define the terrain being investigated need to be defined.

The pressure-sinkage distributions at the wheel–terrain interface are defined by the Preston–Thomas Eq. (4), as it matches well with the data gathered from the penetration tests. The parameters that best fit the measured pressure-sinkage curves (see Fig. 6) were found to be *k* ranging from ≈20 to ≈65 kPa/m and *z*_{m} ranging from 0.3 to 0.54 m. Given this measured variability in the terrain, the parameters of cohesion, angle of internal friction, and shear deformation modulus were tuned to match the average drawbar pull per wheel versus slip and the resistance torque per wheel versus slip given measured data. Lines of drawbar pull versus slip and torque versus slip were generated using the terrain model equations for a single wheel while varying these parameters over the expected range. From this sensitivity study, it was determined that the cohesion parameter in the simulation was to be limited to the lower end of the range defined in the literature for snow [30] to simulate observed immobilizations and match the measured torque. The values for the angle of internal friction, *ϕ*, and the shear deformation modulus, *K*_{s}, were found to have a limited effect on drawbar pull and torque. Therefore, these can be treated as static parameters, with the full range of observed mobility conditions being captured by the measured variability in the pressure-sinkage stiffness parameter, *k*, and with slight changes in cohesion, *c*.

*i*is the slip and

*z*

_{0}is the sinkage that would occur with no excavation and

*S*is the slip sinkage coefficient calculated to be approximately 1.8 [31]. Since this equation was only verified for slips up to 33% [31], it was modified when incorporated into the truth model to limit the additional sinkage at elevated slip ratios, where the added sinkage from slipping, Δ

*z*, can be defined by

*z*

_{0}at 100% slip, rather than ≈180%, as described in Eq. (20). This modified slip sinkage model resulted in an improved approximation of the real data, with the drawbar pull only becoming positive at higher values of

*k*for a narrow band of slip ratios below 50%. At the lower values of

*k*, the sinkage becomes so great that the tractive force that the wheel can apply to the terrain before it fails is not enough to overcome the terrain resistance force. At high values of slip, the slip sinkage becomes so much that even with a relatively high stiffness,

*k*, the drawbar pull still drops below zero. Stepping up the cohesion parameter steps up both the drawbar pull and torque, with its effects more significant at higher slip ratios. Ultimately, the variability of the snow observed in Lebel-sur-Quévillon and the resultant behavior of

*FrostyBoy*can be adequately estimated by varying only

*k*and

*c*. Therefore, the parameters, outlined in Table 2, are used in simulating

*FrostyBoy*operating on low stiffness, poorly bonded, low cohesion snow.

Terrain parameters for simulation | |||||||||
---|---|---|---|---|---|---|---|---|---|

k | z_{m} | c | ϕ | K_{s} | c_{1} | c_{2} | θ_{2} | S | |

Parameter (Units) | $(kPam)$ | (cm) | (kPa) | (deg) | (cm) | – | – | (deg) | – |

Simulation | 20 – 65 | 54 | 0.6 – 1.8 | 9 | 0.03 | 0.18 | 0.32 | −12 | 1.82 |

Terrain parameters for simulation | |||||||||
---|---|---|---|---|---|---|---|---|---|

k | z_{m} | c | ϕ | K_{s} | c_{1} | c_{2} | θ_{2} | S | |

Parameter (Units) | $(kPam)$ | (cm) | (kPa) | (deg) | (cm) | – | – | (deg) | – |

Simulation | 20 – 65 | 54 | 0.6 – 1.8 | 9 | 0.03 | 0.18 | 0.32 | −12 | 1.82 |

#### 3.2.3 Simulations.

Determining the most likely range of terrain parameters by looking at torque and drawbar pull data serve as validation for the performance of each wheel. This could be more rigorously studied using a single-wheel testbed filled with a homogeneous snow with tightly controlled terrain parameters. However, the value of the data obtained during testing is that they were collected in situ, with a four-wheeled robot becoming immobilized on a natural snowpack. This means that validation of the model dynamics of the robot as a whole must be performed by trying to replicate the four-wheeled immobilization signature, rather than solely by the individual wheels operating at a steady state. As discussed in Sec. 3.1, the immobilizations were characterized by two of the wheels, at a diagonal to one another, experiencing a rapid increase in resistance torque until they reach the motor controller current limit, and the other two wheels increasing in rotational speed with a rapid decline of resistance torque until they are almost at a no-load condition. Comparing the torque–slip data points during an immobilization with those of steady-state conditions shows the separation in operating points that occurs and the divergence from the defined torque–slip curves given the range of terrain parameters determined.

Looking at each wheel individually to account for this behavior proved unsuccessful, as no change in the terrain stiffness or cohesion parameters could reproduce the torque–slip operating point during an immobilization, given the sinkages observed. As a result, the dynamic interaction of all four wheels of the rover with the terrain was investigated. Specifically, simulation of transferring the weight of the robot from four wheels to two wheels resulted in the observed effect of high torque and low slip for the weighted wheels and low torque and high slip for the unweighted wheels, with a universal decrease in drawbar pull. This match in wheel behavior resulted in the hypothesis that natural variation in the terrain parameters elicits a shift in weight from all four wheels to only two wheels, ultimately resulting in an immobilization. These two wheels supporting most of the weight must be at a diagonal to one another since this is the only stable configuration. The truth model was then used to investigate how irregularity in the terrain properties could incite this behavior. The simulated terrain parameters below pairs of wheels at a diagonal were arbitrarily changed to the extremes of the ranges outlined in Table 2 (e.g., simultaneous change in stiffness or cohesion of the terrain below the front left and rear right wheels). These unrealistic extreme and sudden changes in terrain properties only resulted in a slight change in operating points for the four wheels and minimal (≈75 N) difference in normal force between the wheels on each axle, with no combination resulting in an immobilization.

*excavation sinkage*term for the wheels that cannot reach their calculated depth due to the rigid chassis. The excavation sinkage is a pseudo-sinkage calculated by subtracting the Δ

*z*term from

*z*

_{0}rather than adding it and is then used to calculate the resistance forces applied to the wheel by the terrain. This has the observed effect of setting up a positive feedback loop between increased slip and weight transfer, as the slipping wheel excavates the terrain that would have supported it. The second parameter added to the model was a wheel sinkage ratio,

*λ*, which is used to calculate the exit angle at the rear of the wheel based on the sinkage at the front of the wheel, instead of keeping this angle constant [32].

*z*in the aforementioned equation,

*θ*

_{2}goes to zero as

*z*goes to zero. The wheel sinkage ratio,

*λ*, was then tuned to ensure that the steady-state torques still matched those observed in the data collected, given the previously defined range of terrain parameters. The exact value of this wheel sinkage ratio in various types of snow needs further study, but for these simulations, a

*λ*value of 0.07 was chosen to give the previously estimated

*θ*

_{2}value of −15 deg at the observed average sinkage of approximately 13 cm.

With these empirically modeled phenomenological effects in place, the rover behavior in the model is consistent with the unique immobilization mode seen in the test data. In addition, the model did not require that the natural, random terrain variation be such that the terrain below two wheels at a diagonal have higher *k* values and the other two wheels simultaneously encounter low cohesion terrain in order to incite a loss of traction. Instead, in the updated model, relatively small variations in the stiffness of the terrain alone can instigate the immobilization behavior consistent with field observations.

## 4 Incipient Immobilization Detection

The challenges in developing a truth model using existing Bekker–Wong and Wong–Reece terramechanics theory and tuning a variety of terrain parameters to match collected data for a wheeled vehicle further reinforce the need for an online mobility assessment algorithm that is independent of these models and parameters. This is proven especially important for operating in a terrain, like snow, with such a wide range of terrain parameters. In snow, a lightweight, wheeled rover must strike a delicate balance between slipping enough to produce the requisite shear displacement to gain traction and slipping so much that slip sinkage results in high terrain resistance, overwhelming the robot’s tractive effort. However, instead of the circuitous method of using estimates of the wheel forces to determine terrain parameters and then calculate a change in operating point, the objective is to look for signatures in the data that can more directly be used to evaluate mobility. The collected data are used to generate a bank of behaviors that are indicative of normal steady-state operation and another set that serve as evidence of incipient immobilization. The MME method is used in real time to decide which hypothesis is the most likely and whether it signifies nominal mobility or immobilizing terrain. In order for this to be valuable for a robot traversing an ice sheet, it must be able to operate in real time. The detection must occur early enough that the rover has time to adapt and maneuver out of the hazardous terrain. Yet, the detection method must be robust enough that it is not triggered by a temporary loss of traction from which the robot can recover. Finally, the bank of models indicative of the two states must be sufficiently universal to predict incipient immobilization even when terrain conditions or robot configurations are different than those directly tested.

### 4.1 Hypothesis Generation.

**H**, to the identity matrix, with the system and input matrices,

**Φ**and

**Γ**, obtained using a prediction error minimization algorithm [34].

This process is easily performed in the System Identification Toolbox in matlab, allowing for rapid testing of input/output combinations with different model orders specified to fit various sections of the time domain data. Each identified system is evaluated based on its goodness of fit with the data from which it is generated using the normalized root-mean-square error cost function. The system parameter values identified from characteristic immobilization segments of data are also compared with those from segments of normal mobility, to assess their uniqueness.

**u**, and robot velocity as the state,

**x**.

**, is made up of a 1×1 matrix,**

*μ***Φ**, and a 1 × 2 matrix,

**Γ**, that modify the state and input, respectively. Hypotheses of nominal mobility are generated from the time series data in which the rover accelerates and reaches steady state. The immobilization hypotheses are generated from the data leading up to the complete loss of forward velocity. Of the 17 immobilization events from the 2 different test sites, only those 6 with the most distinct ramp ups to steady-state operation and subsequent immobilizations were used to generate the hypotheses, with one nominal and one immobilization hypothesis generated for each wheel. The values were plotted for

**Φ**versus

**Γ**(1,1) and

**Γ**(1,2), and the duplicate hypotheses within the same set were removed to limit the number of parallel filters required, leaving 11 nominal hypotheses (

*A*

_{nom}= 11) and 9 immobilization hypotheses (

*A*

_{imm}= 9) for a total of 20 (

*A*

_{total}= 20).

Since the torque–speed curve for a brushless DC motor is well understood and characterized, monitoring the operating torque and wheel speed of the motor and evaluating how this translates to overall rover velocity allows hypothesis selection to run much faster than the 10 Hz rate of data collection without requiring any assumptions about, or solutions to, the equations describing the pressure distribution below the wheels. This method is similar to other condition monitoring techniques, such as the motor current signature analysis (MCSA) used in sensing faults with induction motors [35]. While MCSA uses current spectral analysis to detect anomalies instead of a Bayesian inference algorithm, the basic principle of identifying the operating characteristics of a motor to sense any changes before a rapid decline in performance or a catastrophic failure occurs is embodied by this MME immobilization detection method [36]. It also allows for a robust set of hypotheses to be identified from a small set of immobilization data, meaning that the extensive training data required to train a machine learning algorithm are unnecessary. To prove this, the following sections show how selecting between these 20 hypotheses can be used to reliably detect immobilization not only in the collected data but also in simulated immobilizations generated using the truth model.

### 4.2 Bayesian Multiple Model Estimation Method.

**, is assigned a Kalman filter that determines the likelihood of a measurement given the system parameters within each hypothesis. A posterior probability for each hypothesis is then calculated using Bayes’ rule to combine the likelihood probability of the current time-step with the a priori probability for that hypothesis from the previous time-steps. Typically, each hypothesis is defined by the values of the system matrices (**

*μ***Φ**,

**Γ**, and

**H**) as well as a disturbance input matrix,

**Q**, and a measurement noise matrix,

**R**[37]. For this application, the hypotheses are only defined by matrices

**Φ**and

**Γ**, with

**H**,

**Q**, and

**R**being the same for all hypotheses. Therefore, at each time-step,

*s*, the discrete-time Kalman filter for a given hypothesis,

*μ*_{a}, to calculate an estimate of the state, $x^s$, is given by the sequence:

**r**

_{s}, and the residual covariance matrix,

**S**

_{s}, are calculated with respect to the measured ground speed,

**z**

_{s}, according to

**z**

_{s}, given the estimated state predicted by the hypothesis parameters, $x^s(\mu )$. This distribution is defined by

*A*hypotheses, with the probability of each hypothesis given the measurement defined by

Testing this algorithm on the data collected in Lebel-sur-Quévillon allowed for certain parameters of the Kalman filter to be tuned and for new criteria for the triggering of an immobilization response to be developed. The state estimate depends on the relative weights of the disturbance matrix, **Q**, and the measurement noise matrix, **R**. In general, the ratio of **Q** to **R** defines the relative confidence in the measurements versus the model. Through testing, the best results came when **Q** is set to 20 (m/s)^{2} and **R** is set to 1 (m/s)^{2}. However, both values can vary by $\xb120%$ without any significant impact on detection, showing that the algorithm is not overly sensitive to these parameters. The other tuning mechanism is the adjustment of how frequently a priori probabilities are reset. The recursive Bayesian inference process is designed for convergence, with each new measurement either reinforcing the “belief” that the wheel system is behaving according to the model of a given hypothesis or undermining this “belief.” Gradual convergence would be well suited for slowly changing terrain but is not ideal for the use case for which this incipient immobilization detection method is being designed. The patches of immobilizing terrain are sporadic on the ice sheets of Greenland and Antarctica, and the loss of mobility happens quickly. Therefore, if the *a priori* probabilities have been converging on the hypotheses indicating good mobility for even a few seconds, the sensitivity of the system is reduced. To improve the response time of the system, various methods of flattening the prior distributions at each time-step were tested. Resetting the prior distributions to be uniform naturally resulted in the highest sensitivity and also caused an unacceptable rate of false alarms. Even when requiring that immobilizing hypotheses are chosen for multiple time-steps before a wheel is flagged to be “at risk,” the false alarm rate is still too high, as it would cause the robot to unnecessarily maneuver and expend excess energy to avoid perfectly trafficable snow. Fusing the observed dynamics from all four wheels improves the sensitivity and specificity of predicting mobility with this method.

*A*

_{total}Kalman filters running in parallel, estimating the likelihood of the measured rover velocity given the torque and wheel speed for each hypothesis. These likelihoods are combined with a uniform prior to obtain the respective probability of each hypothesis given the measurement. The average of the probabilities for immobilization hypotheses are subtracted from the average of the probabilities for nominal hypotheses to get a value,

*ξ*, for each wheel, with more negative values signifying higher risk of that wheel becoming immobilized.

*ξ*values from pairs of wheels are summed, requiring two wheels to be “at risk,” before a full immobilization is flagged. The pairs of wheels that most accurately predict an immobilization are those at a diagonal and those on the same axle. Pairing the left-side wheels or the right-side wheels of the robot resulted in false alarms during turning maneuvers. While the robot is highly susceptible to becoming immobilized during a turn, as it needs to overcome the added lateral terrain resistance, a simulated immobilization event during a turn is still detected with these pairings. Therefore, the final detection algorithm calculates $\Xi $ values for each pairing:

*fl*), front right (

*fr*), rear left (

*rl*), or rear right (

*rr*). If any of the four $\Xi $ values drop below a threshold, an immobilization is considered imminent, and the robot needs to take immediate action to maintain mobility.

### 4.3 Detection Results

#### 4.3.1 Collected Data.

The method described earlier was developed and refined using the immobilization event data collected in Lebel-sur-Quévillon. From these data, a threshold value of zero was established for all $\Xi $ values. If any of these four values, in Eqs. (37)–(40), drop below zero, it indicates the likelihoods of the immobilization hypotheses outweigh the likelihoods of the nominal operating hypotheses for that pair of wheels, and the algorithm flags an immobilization. In addition, to prevent false alarms, an immobilization flag is not triggered unless the robot velocity dips below 1 m/s, as the rover’s momentum at higher speeds occasionally allows it to push past immobilizing terrain even if multiple wheels show the signature of an unsustainable motor operating point.

The results of applying the MME method to the data from Lebel-sur-Quévillon are analyzed by looking at two tests. The first test ensures that the algorithm does not trigger in cases where a user initiates a stop and the rover comes to rest due to the motors’ programed deceleration profiles. For the eight events tested in which *FrostyBoy* was commanded to stop rather than becoming immobilized due to the terrain, the algorithm never detected an incipient immobilization. With this simple check accomplished, algorithm’s performance was then assessed on cases where the rover’s velocity did reach a steady-state zero value due to immobilizing terrain. The results of these tests are summarized in Table 3, with Fig. 9 showing an example of when detection occurs (red dotted line) in relation to the rover’s velocity and to each wheel’s rotational speed, torque, and slip. The Supplemental Material on the ASME Digital Collection shows similar plots for all immobilization events observed during the tests near Lebel-sur-Quévillon. For each event, bounds are established for when the rover begins to lose mobility and when the rover has become stuck. While the exact point in time when the robot starts to lose traction is difficult to determine, this is estimated as when the IMU velocity begins to steadily decline to zero and is marked by the first vertical black line in Fig. 9. The second vertical black line marks the time at which the robot’s velocity reaches zero. This allows for algorithm detections that occur outside of these bounds to be excluded from the results’ statistics. Instead, these “early” detections, which do not result in a complete loss of mobility, were analyzed qualitatively, and it was found that the velocity, torque, and wheel speed signatures, with which these “early” flags are associated, are comparable to ones that did result in the robot becoming stuck. In addition, the terrain on which the rover was being tested was already at the edge of trafficability, meaning that if the rover was to encounter this patch of snow on a traverse, it is better to limit how far the rover travels on questionable terrain before raising an alert. So, for all of these cases of “early” detections, it appears that the algorithm functions as intended and does not falsely identify nominal mobility as potentially immobilizing. Also in these “early” cases, the MME method’s final warning was within the established bounds, meaning that the immobilizing terrain that ultimately brought the robot to a halt is identified.

Summary of MME results – Lebel-sur-Quévillon data | |||||
---|---|---|---|---|---|

Event number | Stop distance (m) | Stop time (s) | Velocity at detect (m/s) | Distance to v_{x} = 0 (m) | Time to v_{x} = 0 (s) |

1-1 | 0.54 | 1.2 | 0.97 | 0.45 | 1.1 |

1-2 | 0.88 | 1.5 | 0.86 | 0.38 | 1.0 |

1-3 | 1.34 | 1.8 | 0.45 | Early | Early |

1-4 | 1.58 | 2.6 | 0.91 | Early | Early |

1-5 | 0.46 | 1.0 | 0.82 | 0.15 | 0.7 |

1-6 | 2.66 | 4.7 | 0.58 | 0.67 | 2.1 |

1-7 | 2.31 | 3.7 | 0.95 | 1.26 | 2.8 |

1-8 | 6.07 | 8.0 | 0.98 | 4.42 | 6.6 |

1-9 | 4.69 | 5.2 | 0.99 | Early | Early |

2-1 | 2.64 | 4.0 | 0.88 | 1.51 | 3.0 |

2-2 | 0.38 | 1.3 | 0.84 | Early | Early |

2-3 | 1.95 | 2.4 | 0.65 | Early | Early |

2-4 | 0.78 | 2.4 | 0.56 | 0.61 | 2.1 |

2-5 | 0.67 | 1.7 | 0.70 | Early | Early |

2-6 | 0.45 | 1.1 | 0.70 | 0.21 | 0.8 |

2-7 | 1.00 | 3.4 | 0.29 | 0.28 | 1.8 |

2-8 | 2.29 | 3.6 | 0.82 | 0.99 | 2.2 |

Avg | 1.81 | 2.92 | 0.76 | 0.99 | 2.20 |

Min | 0.38 | 1.00 | 0.29 | 0.15 | 0.70 |

Max | 6.07 | 8.00 | 0.99 | 4.42 | 6.60 |

Summary of MME results – Lebel-sur-Quévillon data | |||||
---|---|---|---|---|---|

Event number | Stop distance (m) | Stop time (s) | Velocity at detect (m/s) | Distance to v_{x} = 0 (m) | Time to v_{x} = 0 (s) |

1-1 | 0.54 | 1.2 | 0.97 | 0.45 | 1.1 |

1-2 | 0.88 | 1.5 | 0.86 | 0.38 | 1.0 |

1-3 | 1.34 | 1.8 | 0.45 | Early | Early |

1-4 | 1.58 | 2.6 | 0.91 | Early | Early |

1-5 | 0.46 | 1.0 | 0.82 | 0.15 | 0.7 |

1-6 | 2.66 | 4.7 | 0.58 | 0.67 | 2.1 |

1-7 | 2.31 | 3.7 | 0.95 | 1.26 | 2.8 |

1-8 | 6.07 | 8.0 | 0.98 | 4.42 | 6.6 |

1-9 | 4.69 | 5.2 | 0.99 | Early | Early |

2-1 | 2.64 | 4.0 | 0.88 | 1.51 | 3.0 |

2-2 | 0.38 | 1.3 | 0.84 | Early | Early |

2-3 | 1.95 | 2.4 | 0.65 | Early | Early |

2-4 | 0.78 | 2.4 | 0.56 | 0.61 | 2.1 |

2-5 | 0.67 | 1.7 | 0.70 | Early | Early |

2-6 | 0.45 | 1.1 | 0.70 | 0.21 | 0.8 |

2-7 | 1.00 | 3.4 | 0.29 | 0.28 | 1.8 |

2-8 | 2.29 | 3.6 | 0.82 | 0.99 | 2.2 |

Avg | 1.81 | 2.92 | 0.76 | 0.99 | 2.20 |

Min | 0.38 | 1.00 | 0.29 | 0.15 | 0.70 |

Max | 6.07 | 8.00 | 0.99 | 4.42 | 6.60 |

Note: The first number of the event number corresponds to the test site. The “early” entries mean that an immobilization was flagged by the algorithm before the start of the immobilization event that actually resulted in the rover’s complete loss of mobility and are not included in the results’ statistics.

For all of the events observed in Lebel-sur-Quévillon, Table 3 shows the wide range of results in how quickly and over what distance the rover can become immobilized. The table also shows that the velocity at the time of detection is generally quite high, with an average value of 0.76 m/s, meaning that a simple threshold on rover velocity would perform significantly worse than the MME algorithm presented. However, the best way to judge the performance of this algorithm is by looking at the distance and time between the immobilization flag and the complete loss of forward velocity for the events without early detection. In these cases, the MME could notify the robot of an impending loss of mobility, on average, ≈1 m, or about 60% of a wheel revolution, before zero velocity. However, this ranged from a mere 0.15 m to a full 4.42 m of warning. While 15 cm is not enough distance to avoid driving into an area of either low cohesion or low stiffness snow, the 0.7 s of response time could drastically reduce the excavation and slip sinkage that occurs at high slips. Also, if the rover has specific control modes to improve its mobility over questionable terrain, or if it is towing a load on a winch, this gives plenty of warning to either alter its operation or pay out cable on the winch to reduce the towed load. The first step is warning the rover of its decreased mobility, and since the MME method triggered an alert for all immobilizations observed with *FrostyBoy* in Lebel-sur-Quévillon, with at least 0.7 s before a complete loss of forward momentum, this method has proven to be a valuable tool in incipient immobilization detection.

#### 4.3.2 Simulated Data.

Due to the unique characteristics of *FrostyBoy*’s immobilizations, from which the system hypotheses were generated, the IDeA, that I present here, is also evaluated by simulation to assess its effectiveness in a different scenario. *FrostyBoy*’s unique immobilizations are attributed to its rigid chassis and the transfer of weight from four wheels to two wheels, as discussed in Sec. 3.2.3. To prove algorithm’s viability for alternate designs of a snow rover, all simulations are run with the assumption that each axle has compliance about the roll axis but otherwise has identical parameters to the *FrostyBoy* rover. The truth model described in Sec. 3.2.1 is used to generate simulated data on which IDeA acts. To incite immobilization, the variability of snowy terrain is modeled by generating random walk models of the terrain stiffness, *k*, and the cohesion, *c*, within the bounds defined by Table 2. Three distinct tests were performed. The first test assessed if the algorithm detects an immobilization for a rover driving straight with all four motors receiving the same open-loop voltage. The second experiment tests the IDeA when all four motors operate under closed-loop speed control. The final test investigates the operation of the rover under closed-loop speed control but with different speed setpoints for the left and right sides, simulating the rover performing a turn.

In all three tests, the IDeA detected the impending immobilization with enough warning that the rover could enter control modes to prevent getting stuck. An example simulation of the first test is shown in the Supplemental Material. In this example, the rover travels at 0.8 m/s when the immobilization is detected, with 0.25 m and 0.6 s until its longitudinal velocity drops to zero. While this is just one example with one random walk model for *k* and *c*, the trend with this test shows that the IDeA detects all simulated immobilizations, but the time and distance between detection and zero velocity decrease compared to the evaluation conducted with field data. This is primarily a result of the rover decelerating much faster in the simulation once entering a high slip condition than that observed in Lebel-sur-Quévillon. This could be an artifact of the simulation, or it could be that in the very different immobilization scenario of all four wheels slipping, the dynamic slip sinkage does cause this rapid increase in terrain resistance as modeled. This requires further investigation in future studies with a newly designed rover.

In operating each motor in the simulation at a constant speed with a closed-loop PI controller and feeding in the same random walk models for *k* and *c* as the aforementioned open-loop cases, the IDeA is more effective at detecting immobilizations earlier and at a longer stopping distance, compared to running in an open-loop mode. For example, in a simulation with closed-loop speed control with identical terrain parameter variation as the first test, the IDeA flags a potential immobilization 0.2 s earlier, meaning that the rover speed at detection is 0.92 m/s and has 0.32 m before its velocity drops to zero. Data from this closed-loop test are found in the Supplemental Material. A comparison of the open-loop versus closed-loop speed control shows this trend across the board, no matter the random walk model that resulted in an immobilization. This is very important since closed-loop speed control is commonly used in autonomous robot control.

The final test adds a turning maneuver to the simulation. During testing in Lebel-sur-Quévillon, the robot was commanded to go straight by giving all four motors the same voltage setpoint in an effort to reduce the variables contributing to its immobilization. In a turn, lateral bulldozing resistance reduces the rover’s ability to gain forward traction and therefore makes it more susceptible to immobilization. Thus, it is most important for the IDeA to warn of a potential loss of mobility during a turn. These simulations are run over a longer timespan, with a different set of *k* and *c* variables than the aforementioned test, to ensure that the robot reached a steady-state yaw rate before an immobilization was initiated. After 4 s of driving straight, the left-side wheels are set to a closed-loop speed setpoint of 3.5 rad/s, while the right-side wheels maintain 5.8 rad/s to turn the rover toward the positive *y*-direction. Of all the simulations run, this study is the furthest departure from the experiments actually performed in the field with *FrostyBoy*. There are no data from Lebel-sur-Quévillon collected during a commanded turn, making it difficult to accurately determine the restoring moment that resists yaw, *M*_{res} from Eq. (??). Also, the truth model is designed for simple longitudinal motion, so there is no explicit *v*_{y} term, and the lateral bulldozing resistance forces are not calculated. Without these factors, the turning maneuver’s negative impacts on mobility are less pronounced in the simulation. However, the model does replicate the negative drawbar pull for the wheels with lower relative speeds (left-side wheels in this case) and the potential for high slips for the wheels with higher relative speeds (right-side wheels). An example of a simulated immobilization that results from this difference between left and right side wheel speed commands can be found in the Supplemental Material. The IDeA is able to predict this loss of mobility in a turn when the simulated rover still has a longitudinal velocity of 0.75 m/s, with 0.8 s and 0.62 m before it comes to a halt. Once again, this proves to be the norm as various random walks of the terrain variables are simulated, as well as different speed setpoints for the right-side and left-side wheels, with the IDeA detecting all immobilizations before the rover loses all momentum.

In these simulations, once a potential loss of traction is detected, the setpoint for all four wheels is set to the same, low-speed value in an attempt to gain full longitudinal traction by commanding the rover to go straight and eliminating the negative effects on mobility resulting from turning. When this simple “traction control” sequence is initiated upon the warning of incipient immobilization during these maneuvers, recovery of mobility is feasible. An example of the effect this had on the rover’s mobility is shown in the Supplemental Material, which uses the same terrain model as the previous turning test, but with both the right- and left-side wheels set to run at 3.5 rad/s once the IDeA warns of an immobilization at 9.3 s. In this simulation, the simple reduction of speed and the command to go straight result in the rover continuing in the longitudinal direction for another 2.7 m beyond where it became immobilized in the previous simulation, where no corrective action was taken (see Fig. 10). The IDeA is also still able to warn of the rover’s final immobilization after recovering from the first. The goal of detailing this example is not to prove out the traction control sequence necessary to navigate low stiffness or poorly bonded terrain when executing a turn. Instead, this example showcases the immobilization detection algorithm’s reliability and its potential in preventing immobilizations, even in scenarios against which it was not specifically ground-truthed.

## 5 Conclusion

The results from applying the IDeA to both real and simulated immobilization data show its effectiveness and reliability in quickly identifying an incipient loss of mobility. The MME only requires torque, wheel speed, and speed over ground measurements as inputs to the model, making it feasible and inexpensive to implement on most rovers. Also, by producing useful results with only 20 hypotheses per wheel and small input matrices, the MME algorithm implementation has a very low computational cost. It is able to run much faster than solving the terramechanics equilibrium equations and faster than similar algorithms involving many more hypothesized terrain parameter sets that each require a separate Kalman filter. Testing this algorithm on field data proves that the dynamics of individual wheels cannot reliably predict impending immobilization, and therefore, the method could not have been developed from single-wheel testbed studies. Instead, detection with a low false alarm rate requires analyzing the response of pairs of wheels and their combined impact on the velocity of the rover. Finally, testing this algorithm on simulations shows that this method is not only applicable to the data set from which it was generated. Instead, the promising results from flagging immobilizations generated by the truth model in scenarios outside the scope of the experiments with *FrostyBoy* in Lebel-sur-Quévillon indicate that this method could be applied to rovers that are of a different design, that encounter more variation in the terrain, and that operate under a different control scheme. The truth model, the detection algorithm, and the data from the immobilization events are made available on GitHub [38] to promote further research of this technique. Additional studies need to be performed to test this detection method with the configurations and robot geometries that have only been simulated. In addition, for a wider set of nominal and immobilization hypotheses to be generated from a small data set in the future, the truth model should be further refined and validated from *in situ* testing. Ultimately, the results from this study show that the IDeA is a simple, robust, and reliable method for detecting incipient immobilization of a rover operating in snow.

## Footnote

## Acknowledgment

This work was supported by NSF’s program in the division of Civil, Mechanical and Manufacturing Innovation (CMMI) for Dynamics, Control and System Diagnostics (DCSD) (Award No. 1824687).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The data and information that support the findings of this article are freely available.^{2}

## Nomenclature

*a*=hypothesis number

*b*=wheel width

*c*=cohesion

*g*=acceleration due to gravity

*h*=height of center of mass above axles

*i*=slip ratio

*j*=shear displacement

*k*=terrain stiffness

*m*=mass of rover

*n*=exponent, modifying sinkage in pressure-sinkage equations

*p*=pressure

*q*=distance step number for random walk models

*s*=time-step number in discrete state space model

*t*=time

*z*=sinkage, distance from terrain surface to BDC of wheel

**r**=measurement residual vector

**u**=input vector

**x**=state vector

**y**=output vector

**z**=measurement vector

*S*=slip sinkage coefficient

*T*=torque applied by wheel

*W*=weight on each wheel

**H**=state space output matrix

**K**=Kalman filter gain

**Q**=process disturbance covariance matrix

**R**=measurement noise covariance matrix

**S**=residual covariance matrix

- $h^$ =
Grouser height normalized to wheel radius

- $z^$ =
sinkage normalized to wheel radius

- $v\u02d9x$ =
acceleration along longitudinal axis

- $x^\u2212$ =
a priori state estimate

- $x^+$ =
posteriori state estimate

*c*_{1}=coefficient 1 used to determine location of

*θ*_{m}*c*_{2}=coefficient 2 used to determine location of

*θ*_{m}*d*_{f}=distance from center of mass to front axle

*d*_{r}=distance from center of mass to rear axle

*k*_{c}=cohesive modulus, Bekker pressure-sinkage model

*k*_{ϕ}=friction modulus, Bekker pressure-sinkage model

*r*_{w}=wheel radius

*t*_{w}=track width

*v*_{L}=velocity of left wheels along longitudinal axis

*v*_{R}=velocity of right wheels along longitudinal axis

*v*_{x}=velocity along longitudinal axis

*z*_{0}=sinkage, without taking into account slip sinkage effect

*z*_{actual}=sinkage, taking into account slip sinkage effect

*z*_{f}=sinkage of front axle wheels

*z*_{m}=maximum sinkage, Preston–Thomas pressure-sinkage model

*z*_{r}=sinkage of rear axle wheels

*z*_{t}=depth of transition point of forward region and rear region flow regimes

*A*_{imm}=number of immobilization hypotheses

*A*_{nom}=number of nominal mobility hypotheses

*A*_{total}=sum of number of immobilization and nominal mobility hypotheses

*B*_{ω}=mechanical damping term

*F*_{x}=drawbar pull

*F*_{xfl}=drawbar pull of front left wheel

*F*_{xfr}=drawbar pull of front right wheel

*F*_{xrl}=drawbar pull of rear left wheel

*F*_{xrr}=drawbar pull of rear right wheel

*F*_{z}=normal force on wheel

*F*_{zf}=total normal force on both wheels of front axle

*F*_{zfl}=normal force on front left wheel

*F*_{zfr}=normal force on front right wheel

*F*_{zr}=total normal force on both wheels of rear axle

*F*_{zrl}=normal force on rear left wheel

*F*_{zrr}=normal force on rear right wheel

*I*_{zz}=moment of inertia of rover about the

*z*-axis at its center of mass*I*_{ω}=wheel moment of inertia about its center of rotation

*K*_{s}=shear deformation modulus, poorly bonded shear stress–shear displacement model

*M*_{res}=restoring moment

*T*_{fl}=torque applied by front left wheel

*T*_{fr}=torque applied by front right wheel

*T*_{r}=resistive torque

*T*_{rl}=torque applied by rear left wheel

*T*_{rr}=torque applied by rear right wheel

*T*_{rfl}=resistive torque applied by terrain to front left wheel

*T*_{rfr}=resistive torque applied by terrain to front right wheel

*T*_{rrl}=resistive torque applied by terrain to rear left wheel

*T*_{rrr}=resistive torque applied by terrain to rear right wheel

- $Ps\u2212$ =
a priori state estimate covariance matrix

- $Ps+$ =
posteriori state estimate covariance matrix

*β*=pitch

*γ*=specific weight of snow or soil

**Γ**=state space input matrix

- Δ
*z*= additional sinkage due to slip

*θ*=angle from BDC of wheel

*θ*_{1}=angle from BDC of wheel to surface of the terrain in front region

*θ*_{2}=angle from BDC of wheel to where wheel exits the terrain in rear region

*θ*_{m}=angle from BDC of wheel to the point of max pressure

*λ*=wheel sinkage ratio

**Λ**=state space disturbance matrix

=*μ*hypothesis matrix

*ξ*=wheel nominal mobility value

- $\Xi $ =
wheel pairing nominal mobility value

*ρ*=roll

*σ*=normal stress

*τ*=shear stress

*ϕ*=angle of internal friction

**Φ**=state space system matrix

*ψ*=yaw

- $\psi \u02d9$ =
yaw rate

- $\psi \xa8$ =
yaw angular acceleration

*ω*=angular wheel speed

*ω*_{fl}=angular wheel speed of front left wheel

*ω*_{fr}=angular wheel speed of front right wheel

*ω*_{rl}=angular wheel speed of rear left wheel

*ω*_{rr}=angular wheel speed of rear right wheel

- $\omega \u02d9fl$ =
angular wheel acceleration of front left wheel

- $\omega \u02d9fr$ =
angular wheel acceleration of front right wheel

- $\omega \u02d9rl$ =
angular wheel acceleration of rear left wheel

- $\omega \u02d9rr$ =
angular wheel acceleration of rear right wheel