2.1 Uncertainty Modeling

Robot motion in real-world environments is inherently uncertain because environmental disturbances disrupt the mapping between commanded inputs and physical movement. In differential-drive systems, even small systematic errors such as unequal wheel diameters or imprecise wheelbase measurements accumulate into drift over time and degrade the reliability of open-loop control [19]. Cardinaels et al. [14] address environmental uncertainty in predictive teleoperation through a worst-case deterministic formulation. Their approach represents uncertainty by constructing a bounded envelope over predicted motion under asymmetric disturbance assumptions, where one wheel is subjected to worst-case conditions while the other operates normally. This yields a conservative geometric bound on lateral deviation arising from accumulated disturbances.

Beyond deterministic bounds, probabilistic modeling has become central to characterizing the stochastic nature of motion disturbances. Thrun et al. [20] formalized this through velocity and odometry motion models in which Gaussian noise is injected into translational and rotational components of commanded motion. This produces a distribution over successor poses rather than a single predicted state and established the foundation for treating motion commands as stochastic processes.

Accurately propagating these distributions through nonlinear motion models, however, introduces significant computational complexity that limits their direct use in real-time systems. The Unscented Transform (UT) introduced by Julier and Uhlmann [21] addresses this challenge by approximating the propagation of mean and covariance through nonlinear transformations using a deterministic set of carefully chosen sigma points. Unlike linearization-based approaches, UT captures higher-order effects of nonlinearity without requiring explicit Jacobians, making it particularly suitable for real-time state estimation in robotics.

Building on these principles, data-driven approaches have been developed to capture terrain-induced disturbances in off-road and planetary contexts. Trivedi et al. [22] employed Gaussian Process Regression augmented with sigma-point transforms to model nonlinear tire–terrain interactions, generating full predictive distributions over future poses. Cunningham et al. [23] used Gaussian processes to estimate wheel slip for planetary robots, with uncertainty bounds that tighten as local terrain observations accumulate. Probabilistic state estimation under communication delay has also been explored. Das and Dobie [24] modeled the distributional uncertainty of delayed sensor measurements within an augmented-state Kalman filter, allowing the estimator to reason explicitly about temporal misalignment.

While these models improve predictive fidelity, the central challenge in delayed teleoperation is not solely the accuracy of the underlying uncertainty representation. Instead, it concerns whether the resulting predictions and their associated uncertainty can be externalized in a form that operators can perceive, interpret, and act upon within the tight temporal constraints of real-time control.

2.2 Visual Encoding

Probabilistic models of robot motion provide mathematically precise descriptions of uncertainty, but making such information actionable for human operators requires careful interface design [25]. Probability estimates are often difficult to interpret in real time, particularly under time pressure or for non‑expert users [26, 27]. Visualization therefore serves as a critical translation layer, converting probabilistic information into perceptual forms that can be integrated directly with an operator’s understanding of the robot and its environment.

Uncertainty visualization is a broad and heterogeneous design space rather than a single representational solution [28, 29]. Foundational surveys, such as Pang et al. [30], catalog a wide range of techniques, including fuzzy boundaries, overlays, and variations in color, opacity, texture, and geometry, highlighting that uncertainty can be encoded in fundamentally different ways. A central distinction in this space lies between bounded encodings, which communicate only the range of possible outcomes, and distributional encodings, which convey how probability mass is distributed within that range.

Bounded encodings, formalized by Olston and Mackinlay [31], represent uncertainty using explicit upper and lower limits without indicating the likelihood of intermediate outcomes. Predictive visualizations commonly used in teleoperation, including the worst-case envelope described above, fall into this category. These representations are visually simple and computationally efficient. However, they do not communicate how likelihood varies within the boundary, potentially obscuring structure that is important for anticipating likely outcomes.

Distributional encodings aim to make probability distributions perceptually apparent. Techniques include gradient-based fills, density and violin plots, contour lines, and dynamic ensemble-based cues [30, 32, 33, 29]. Empirical studies show that such representations communicate information that bounded displays cannot. Seipp et al. [34] found that opacity-based encodings support more accurate interpretation of areal uncertainty than alternatives such as blur, grids, texture, or contour lines, particularly for non-expert users. Complementary work on Value‑Suppressing Uncertainty Palettes shows that reducing color resolution in regions of high uncertainty encourages more cautious decision‑making [35].

Colley et al. [36] operationalized these principles in the context of autonomous vehicles, implementing a cone-shaped trajectory visualization that widens to encode growing positional uncertainty, with color and transparency indicating likelihood. Their results confirmed that color is particularly effective at conveying distributional uncertainty. Taken together, these findings underscore that visual encoding fundamentally shapes whether probabilistic information becomes actionable insight, motivating the shift from bounded to distributional representations in predictive display design.

2.3 Human Interpretation

Whether a visualization of uncertainty supports effective control depends not only on its design, but on the cognitive processes operators bring to bear when reading and acting on it under delay. When immediate feedback is unavailable, operators must maintain internal mental models to anticipate the robot’s unobserved state [6, 12, 4, 5]. Constructing and updating these models increases mental workload, as reflected in studies reporting elevated cognitive demand, heightened stress, and conservative behavioral strategies such as move‑and‑wait under delayed conditions [6, 5, 7, 8]. Predictive visualizations aim to externalize uncertainty and reduce this burden, but their effectiveness depends critically on how operators perceive and interpret probabilistic information [37].

Boukhelifa and Duke [38] argue that uncertainty visualizations often fail in practice due to perceptual and cognitive confusion introduced by the depiction itself. They identify the aggregation of heterogeneous uncertainty sources as a particular challenge, noting that combining confidence information from multiple origins can introduce a secondary layer of uncertainty that undermines the reliability of the display. Matzen et al. [37] reinforce this point empirically, showing that seemingly minor differences in how uncertainty is depicted can lead to substantially different patterns of risk perception and decision‑making, underscoring that depiction choices are not neutral.

Empirical work in related domains further demonstrates that visual encodings of likelihood shape decision‑making, risk perception, and trust. Colley et al. [36, 17] show that visualizing uncertainty can improve users’ anticipation of likely outcomes and calibrate trust in the underlying system, while Bancihon et al. [39] find that the choice of uncertainty encoding modulates both perceived risk and risk‑taking behavior. Large‑scale studies by Bles et al. [40] indicate that communicating uncertainty bounds does not significantly undermine trust in the information source, and Reyes et al. [41] extend this finding to AI‑based decision support, showing that uncertainty visualization can meaningfully increase trust among initially skeptical users. Although these studies do not focus exclusively on teleoperation, they suggest a transferable design principle: making likelihood information perceptually salient can help users distinguish likely outcomes from merely possible ones. In the present work, this principle motivates the use of layered probability bands, while its operational effect in teleoperation remains a question for future user evaluation.

Taken together, prior work highlights the need to move beyond bounded representations toward predictive displays that convey the full distribution of possible outcomes. This shift emphasizes not only probabilistic modeling, but also how uncertainty is encoded and interpreted by human operators. Distribution-aware visualizations have the potential to reduce cognitive load, improve decision-making under delay, and avoid overly conservative control strategies.

3.1 Robot Motion

UNITE simulates a TurtleBot3 Waffle Pi, a differential‑drive mobile robot whose motion follows the standard kinematic model. The operator provides discrete movement commands (forward, backward, left, right), which are translated by the simulator into left and right commanded wheel velocities, denoted \(v_{L}\) and \(v_R\). At each simulation step \(dt\), the simulator updates the robot’s pose and wheel velocities. For implementation purposes, UNITE stores the robot’s position \((x,y,z)\), heading \(\theta\), wheel velocities \((v_L,v_R)\), timestep \(dt\), and current input. The probabilistic model developed in this work primarily uses the planar pose \((x,z,\theta)\) and the wheel velocities \((v_L,v_R)\). The robot’s motion is determined by its linear and angular velocities. The linear velocity \(v\) is computed as the average of the wheel velocities: \[v = (v_L + v_R) / 2\] and the angular velocity \(\omega\) is computed from their difference, scaled by the robot’s wheelbase distance \(b\): \[\omega = (v_R - v_L) / b\] These velocities are then integrated over the timestep \(dt\) to update the robot’s pose: \[ x_{t+1} = x_t + v \cos(\theta_t)\, dt, \] \[ z_{t+1} = z_t - v \sin(\theta_t)\, dt, \] \[ \theta_{t+1} = \theta_t + \omega\, dt. \]

3.2 Communication Delay

To reproduce the temporal separation inherent to lunar teleoperation, UNITE incorporates a fixed one-way communication delay of 2.56 s. This value corresponds to the minimum signal propagation time between Earth and the Moon [21], thereby grounding the simulation in a realistic operational context.

In practice, the delay is implemented using a time‑stamped command queue. Each operator input is stored in this queue together with the simulation time at which it is scheduled to be executed. Rather than being applied immediately, commands are deferred until the full 2.56 s delay has elapsed.

During this waiting period, the queued inputs are also used to compute the predicted lookahead window. This ensures that the robot’s anticipated future trajectory is always derived from the same delayed command stream that governs its actual execution.

3.3 Environmental Disturbances

Beyond communication delay, UNITE incorporates a diverse set of environmental disturbances that reflect the variability encountered in planetary and lunar environments. These disturbances affect the robot’s wheel velocities before each kinematic update, causing deviations between commanded and executed motion. They include:

1. Wheel Slip: Reduced traction on rough or loose surfaces lowers the effective wheel speeds, especially at higher velocities.
2. Motor Variation: Time-dependent, deterministic variations in motor response independently modulate the left and right wheel velocities, introducing small but continuous asymmetries.
3. Wheel Radius Variation: Slight differences in wheel radius introduce a slow, oscillatory turning tendency even under symmetric commands.
4. Terrain Vibration: Spatially smooth terrain irregularities generate fluctuations in wheel motion.
5. Encoder Noise: Additional variability is applied to wheel velocities to reflect imperfections.
6. Slope Effects: Inclined terrain alters effective wheel velocities, causing downhill drift or resistance when moving uphill.

Disturbances (1) and (6) are directly governed by the terrain at the robot’s position. Wheel slip scales with both local surface roughness and the robot’s velocity, whereas slope effects are determined by the terrain gradient and the robot’s heading, introducing direction-dependent drift or resistance. Disturbances (2) and (3) are deterministic and time-dependent: motor variation is modeled through sinusoidal functions, while wheel-radius variation is represented by an anti-correlated sinusoidal signal applied to the left and right wheels. Disturbances (4) and (5) are derived from Perlin noise fields sampled at the robot’s location, generating spatially coherent vibration and encoder fluctuations that evolve smoothly as the robot traverses the environment. Although the resulting motion may appear irregular, all disturbance sources within UNITE’s physics model are deterministic, meaning that identical inputs consistently produce identical outcomes. This contrasts with stochastic uncertainty models, where disturbances are treated as random variables rather than fixed, repeatable dynamics.

To account for these environmental influences, the nominal wheel velocities \(v_{L}\) and \(v_R\), computed directly from operator input, are subsequently processed by the physics model. This yields modified wheel velocities that capture the combined effects of control input and terrain-dependent disturbances.

4.1 Modeling of Environmental Disturbances

To compute the mean trajectory and the corresponding sigma-bands, we separate disturbances in the physics model into two categories based on their characteristics.

Fully deterministic. As described in Section 3.3, we treat slip, motor variation, wheel radius variation, and slope bias as deterministic effects. These are modeled using functions whose future values are known exactly at the time of prediction. For instance, motor and wheel radius variations are represented as fixed-frequency sine waves that can be evaluated directly at any future timestep. The slope bias depends only on the terrain gradient and the rover’s heading, both of which are known in advance.

Slip is also treated as deterministic in this formulation. Given the current speed and terrain roughness, it produces a symmetric reduction in wheel speed that can be computed exactly. Since all of these effects can be determined for future timesteps without uncertainty, they are not treated as stochastic terms. Instead, they are directly incorporated into the mean wheel velocity.

Temporally correlated. The remaining disturbances, namely encoder noise and vibration noise, are treated as stochastic. As discussed in Section 3.3, they are modeled using Perlin noise fields, meaning their exact future values cannot be predicted. Unlike white noise, Perlin noise is smooth over time, so consecutive values are correlated rather than independent.

To model this temporal correlation in a tractable way, we use a first-order Gauss-Markov process. This is a natural choice because it captures temporal dependence through an exponentially decaying influence of past values.

4.2 Forward Propagation

At each update step, the mean trajectory and the sigma-bands are computed together in a single forward pass through the sequence of delayed control inputs. First, deterministic disturbances are applied to the commanded velocities to obtain the mean wheel velocities $\bar{v}_{L,k}$ and $\bar{v}_{R,k}$. These values are then propagated through the differential-drive kinematics to update the mean pose.

The sigma-bands require tracking how uncertainty evolves over the lookahead horizon. Because the differential-drive dynamics are nonlinear, linearization-based methods such as the Extended Kalman Filter can introduce errors that accumulate over time. Instead, the UT represents uncertainty using a small set of sample points. These points are propagated through the exact nonlinear motion model, and the resulting outputs are used to reconstruct the updated mean and covariance. This approach avoids linearization errors and provides a more accurate estimate of how uncertainty evolves over time, while still remaining computationally efficient.

Since vibration and encoder noise cannot be predicted exactly, their uncertainty must be propagated along with the kinematic state. For this reason, the four Gauss-Markov noise states are included directly in the state vector, resulting in a seven-dimensional representation: \[ \mathbf{x}_k = \bigl(x, z, \theta, n_{\text{vib},L}, n_{\text{vib},R}, n_{\text{enc},L}, n_{\text{enc},R}\bigr)_k. \]

This formulation allows the UT to propagate both pose and noise uncertainty in a single step, ensuring that the accumulated effect of temporally correlated disturbances is reflected in the sigma-bands as the prediction horizon increases.

We maintain a corresponding $7 \times 7$ covariance matrix $\Sigma_k$. At each step, it is propagated using the UT by first computing its Cholesky decomposition, which yields a lower-triangular matrix $L_k$ such that \[ \Sigma_k = L_k L_k^\top, \]

The columns of $L_k$ are then used to construct $2n = 14$ symmetric sigma points: \[ \mathcal{X}^{(i)}_k = \bar{\mathbf{x}}_k \pm \sqrt{n}\,(L_k)_i, \quad i = 1, \ldots, n, \] where $n = 7$ is the state dimension and $(L_k)_i$ denotes the $i$-th column of $L_k$.

The standard UT typically includes an additional sigma point at the mean, resulting in $2n+1$ points. In this work, we use a simplified parameterization with $\alpha = 1$ and $\kappa = 0$, which sets the scaling parameter $\lambda = \alpha^2(n+\kappa) - n$ to zero. With this choice, the central point has zero weight and is omitted, leaving $2n$ equally weighted points with $W_i = \tfrac{1}{2n}$.

Each sigma point is then propagated through the nonlinear differential-drive kinematics, and the updated covariance is reconstructed as \[ \Sigma_{k+1} = \sum_{i=1}^{2n} W_i \left(\mathcal{X}^{(i)}_{k+1} - \bar{\mathbf{x}}_{k+1}\right) \left(\mathcal{X}^{(i)}_{k+1} - \bar{\mathbf{x}}_{k+1}\right)^\top + Q_k, \]

The matrix $Q_k$ represents process noise added at each propagation step. It has two components. First, noise is injected into the vibration and encoder states to prevent their uncertainty from artificially shrinking over long horizons. Second, a residual pose-level noise term is included to capture effects not modeled in the deterministic dynamics, such as unmodeled slip and small terrain irregularities. This pose-level term is projected along the rover’s heading and wheelbase, preserving correlations between $x$, $z$, and $\theta$.

After obtaining \(\Sigma_{k+1}\), the lateral uncertainty \(\sigma_{\text{lat}}\) is computed by projecting the positional covariance onto the lateral direction vector \((\sin\theta_k,\ \cos\theta_k)^\top\): \[ \sigma_{\text{lat}} = \sqrt{\sin^2\theta_k\,\Sigma_{xx} + 2\sin\theta_k\cos\theta_k\,\Sigma_{xz} + \cos^2\theta_k\,\Sigma_{zz}}, \]

Finally, the \(\pm 1\sigma\), \(\pm 2\sigma\) and \(\pm 3\sigma\) bounds are constructed symmetrically around the mean trajectory along this lateral direction, forming the Probabilistic Envelope visualized during teleoperation.

6.1 Computational Efficiency

The main motivation for adopting the Unscented Transform (UT) over Monte Carlo (MC) rollouts was real-time feasibility. MC was still included because it provides a straightforward way to model non-Gaussian uncertainty without relying on Gaussian assumptions. To evaluate both approaches, we benchmarked the execution time required to compute the lookahead envelope on a per-frame basis.

In the MC approach, the algorithm iterates over all inputs stored in the delayed command queue. For each input, a set of simulated trajectories is propagated forward in time. All samples start from the same initial state but evolve differently due to stochastic perturbations, which are decorrelated using unique Perlin-noise offsets per sample. This allows the method to avoid Gaussian assumptions and produce a trajectory distribution that better reflects the variability of the underlying dynamics.

At each prediction step, the lateral deviation of each sample is computed by projecting its displacement onto the lateral axis defined by the mean heading. These deviations are then sorted to form an empirical distribution, from which percentile-based bounds are extracted. The inner band corresponds to the 25th–75th percentiles (50% of samples), the middle band to the 5th–95th percentiles (90% of samples), and the outer band to the 1st–99th percentiles (98% of samples). These percentile offsets are then mapped back into spatial coordinates to construct the final envelope boundaries.

Both methods were evaluated under controlled conditions. The simulator generated predictions over 500 frames, each using a horizon of 128 input steps determined by the communication delay and the physics update rate. For the MC method, the computation was repeated with different numbers of samples to study how runtime scales with trajectory count.

The experiments were conducted under three physics profiles representing a lunar-like rough surface, a smooth flat surface, and loose sand. Each profile captures a different operating regime relevant to off-road rover dynamics. The parameter configurations are listed in Table 1.

Figures 2 and 3 compare elapsed time across physics profiles and methods. The red dashed line indicates the 8.33 ms per-frame budget at 120 FPS, which must accommodate rendering, physics, input handling, and envelope computation. Exceeding this budget by the envelope alone makes frame drops unavoidable, while staying below it does not guarantee real-time performance due to other Unity overhead.

The UT method averages 800 µs per frame across all physics profiles, whereas MC ranges from 10 ms (100 samples) to 63 ms (550 samples), exceeding the frame budget in all configurations. UT remains below the threshold in nearly all cases, with only three outliers slightly above it.

Figure 4 shows that elapsed time scales linearly with sample count across all physics profiles. The Moon profile incurs a consistently higher cost than FlatSurface and LooseSand, a difference that grows from approximately 1.5 ms at 100 samples to 10 ms at 550 samples.

6.2 Predictive Calibration and Containment

To verify that the predicted uncertainty envelopes accurately reflect the physical disturbances present in the simulator, we evaluate the containment rate of the Probabilistic Envelope. If the probabilistic model is well calibrated, the robot’s actual trajectory should fall within the $1\sigma$ boundary approximately 68% of the time, within the $2\sigma$ boundary approximately 95%, and within the $3\sigma$ boundary approximately 99.7% of the time, consistent with the expected containment probabilities of a Gaussian distribution.

To evaluate containment, three quantities were logged at each simulation frame. First, the robot’s actual position in world space. Second, the predicted mean trajectory produced at that frame, storing the expected $(x, z)$ position for each step in the prediction horizon alongside the originating frame index and the corresponding target frame index. Finally, the left and right boundary points of the $1\sigma$, $2\sigma$, and $3\sigma$ sigma-bands were recorded for each horizon step as $(x, z)$ pairs, each associated with its from-frame and target-frame indices.

The experiment used the moon physics preset designed to reproduce key disturbance characteristics of low-gravity, irregular terrain. The corresponding physics parameters are listed in Table 2.

The robot was driven manually along a single constrained path, since the terrain geometry allowed no alternative routes. The trajectory covered 15.05 m over 167.8 s.

The predictive behavior of the mean trajectory is illustrated in Figure 5(a). At each frame, the model generates a full forward prediction from the robot’s current state independently, without using any temporal memory of previous predictions. To keep the visualization readable, 100 frames are uniformly sampled from the recorded sequence.

For each sampled frame, the predicted mean positions over the horizon are connected to form a short trajectory segment starting from the robot’s actual position at that frame. When placed along the full trajectory, these segments show how the model’s predictions evolve as new observations become available. Segments that closely follow the actual trajectory indicate accurate mean predictions, while segments that deviate highlight regions where prediction reliability decreases. As shown in Figure 5(a), the predicted segments remain closely aligned with the robot’s trajectory throughout the run.

While the mean trajectory captures expected motion, it does not provide information about uncertainty. Figure 5(b) therefore evaluates whether the robot’s actual position lies within the $3\sigma$ bounds of the Probabilistic Envelope at each frame. For each prediction, a lateral containment test is performed: the robot’s actual position is projected onto the lateral axis defined by the predicted mean and the left $1\sigma$ boundary, producing a signed lateral displacement relative to the predicted mean.

A frame is marked green if at least 95% of the predictions reaching that frame contain the robot within the corresponding $3\sigma$ lateral interval, and red otherwise. As shown in Figure 5(b), the trajectory is mostly green, with a short red section near the end of the run.

To quantify these results, Table 3 reports containment rates for the $1\sigma$, $2\sigma$, and $3\sigma$ bands. Since a new prediction is generated at every frame and each prediction spans multiple future horizon steps, the evaluation produces a large number of prediction-observation pairs obtained by merging actual positions, mean predictions, and sigma envelopes across matching frame–horizon indices. For each pair, the robot’s actual future position is expressed relative to the predicted mean and projected onto a locally defined lateral axis derived from the $1\sigma$ boundary geometry. The containment test is then performed by computing the signed lateral offset of the actual position with respect to the mean and checking whether it lies between the projected left and right sigma bounds for each confidence level. The containment rate is computed as the fraction of pairs where this offset lies within the corresponding $1\sigma$, $2\sigma$, and $3\sigma$ intervals.

The $1\sigma$ and $2\sigma$ containment rates exceed their theoretical Gaussian expectations, with observed rates of 80.41% against an expected 68.27% for $1\sigma$, and 97.69% against an expected 95.45% for $2\sigma$. The $3\sigma$ containment rate of 99.28% falls marginally below the theoretical value of 99.73%.