See the Future, Classify the Present: Stochastic Latent World Modeling for Hyperspectral Images

Introduction

$\text{H}$yperspectral image (HSI) classification (HSIC) aims to assign semantic labels to pixels by leveraging high-dimensional spectral signatures together with spatial context [1]. The rich spectral resolution of hyperspectral imagery enables fine-grained material discrimination [2]; however, it also introduces significant challenges, including high dimensionality, strong inter-band correlation, and severe scarcity of labeled samples [3]. Despite substantial advances, most existing HSIC approaches remain fundamentally based on static spatial-spectral representations [4], where classification is performed from a single-stage encoding of local patches [5]. This paradigm limits the ability of models to capture higher-order dependencies and structured relationships beyond immediate context, particularly in low-label regimes [6].

Limitations of Existing Approaches. Traditional HSIC methods, including hand-crafted feature pipelines and early deep models, represent each patch as a fixed embedding that is directly mapped to a class label [7]. Convolutional neural networks (CNNs) and spatial-spectral hybrid models improve representation capacity through hierarchical feature extraction and local context aggregation [8], [9]. However, these models remain inherently feed-forward and static: the learned representation is a single point estimate in latent space [10]. Such formulations do not explicitly encode relationships between latent states, nor do they model how representations could evolve under structured transformations [11]. Consequently, they fail to exploit the continuity and smoothness inherent in hyperspectral data [12], [13].

Attention-Based Context Modeling Without Structure. Transformer-based models have been introduced to capture long-range dependencies across spectral bands and spatial locations [14, 15, 16]. While multi-head self-attention (MHSA) improves global context aggregation [17], it operates as a weighted averaging mechanism over tokens within a single forward pass [18]. The resulting representation is still a one-step aggregation without explicit constraints on latent organization [19, 20, 21]. In particular, attention does not impose transition structure or predictive consistency across representations, which limits its ability to enforce physically meaningful relationships such as spectral continuity and spatial coherence [22, 23, 24].

Self-Supervised Learning Without Latent Transitions. Self-supervised approaches, including masked reconstruction [7], [25], contrastive learning [26], [27], and generative modeling [28], [29], aim to reduce dependence on labeled data by introducing auxiliary objectives [30]. These methods encourage representations to preserve information or maintain invariances [31]. They do not explicitly model dependencies between latent states, nor do they enforce how representations should transform under structured variations. As a result, they provide regularization but not structured predictive modeling [32].

World Models and Their Limitations. Recent advances in world models and latent dynamics learning [33], particularly in reinforcement learning and sequential modeling, aim to learn compact representations that evolve through predictive transitions [34, 35, 36]. These methods typically rely on explicit temporal observations and learn state transition dynamics

M. Ahmad and A. Sohaib are with the Khanpur Institute of Technology, Khanpur 64100, Pakistan. (e-mail: [email protected]; [email protected]).

G. Vivone is with the National Research Council - Institute of Methodologies for Environmental Analysis (CNR-IMAA), 85050 Tito, Italy (e-mail: [email protected]).

R. Hänsch is with the Department of SAR Technology, German Aerospace Center (DLR), Germany. (e-mail: [email protected]).

M. Mazzara is with the Institute of Software Development and Engineering, Innopolis University, Innopolis, 420500, Russia. (e-mail: [email protected]).

D. Hong is with the School of Automation, Southeast University, Nanjing, China. (e-mail: [email protected]) conditioned on actions or sequential inputs [37, 38]. While effective in domains with well-defined temporal structure, such formulations are not directly applicable to hyperspectral imagery, which is inherently static and lacks explicit temporal supervision [39].

Moreover, existing world models assume that transitions correspond to physically meaningful state changes observed in time [40, 41, 42, 43, 44], whereas hyperspectral data primarily exhibit structural correlations rather than temporal evolution. Applying standard world modeling frameworks to HSIC would therefore require artificial sequence construction or surrogate temporal assumptions, which may introduce bias and instability. In addition, world models focus on generative prediction in observation space, rather than enforcing structured organization within latent representations for discriminative tasks [45, 46, 47]. Consequently, they do not directly address the need for latent predictive regularization in static spatial-spectral data [48, 49].

Absence of Latent Predictive Structure. A key limitation across these categories is the lack of latent predictive modeling tailored to HSIC. Hyperspectral data exhibit strong intrinsic structure, including smooth spectral variation across bands and spatial consistency across neighboring pixels. However, current models do not capture these properties as transitions in latent space. Without modeling such transitions, representations remain locally optimized but globally unstructured. This restricts the ability to generalize under limited supervision, as the model cannot reason beyond the immediate encoded observation.

Proposed Perspective

1) Proposed Perspective: To overcome these limitations, this work introduces stochastic latent world modeling (LWM), a trajectory-aware representation learning framework for HSIC. The central idea is to move beyond static encoding and instead learn representations that are structured through latent evolution. Rather than treating the representation of a patch as a fixed embedding, LWM models how this representation can transform under a learned stochastic process in latent space.

Unlike classical world models that rely on explicit temporal observations, the proposed formulation operates entirely in representation space. A latent environment state is first constructed from fused spatial-spectral features. A stochastic transition model then generates a sequence of latent states via multi-step rollouts, parameterized as Gaussian transitions with recurrent refinement. These rollouts represent plausible variations of the underlying representation, acting as a form of structured regularization rather than modeling physical time.

To extract decision-relevant information, a rollout planning mechanism adaptively aggregates latent states using learned importance weights. This produces a trajectory-aware representation that complements the base transformer encoding. The final prediction is obtained by combining the static and trajectory-based representations, enabling the model to leverage both immediate evidence and structured latent evolution.

The framework is trained with complementary objectives that jointly enforce reconstruction fidelity, latent dynamics consistency, stochastic regularization, and physically meaningful constraints. In particular, spectral and spatial dynamics losses explicitly model structured relationships within latent space, while stochastic rollout regularization enforces probabilistic coherence across predicted transitions. This unified formulation enables the model to learn representations that are not only discriminative but also structured and robust. The main contributions of this work are summarized as follows:

• Trajectory-Aware Representation Learning for HSIC: We reformulate HSIC as a trajectory-aware latent representation problem, moving beyond static spatial-spectral encoding toward structured latent reasoning.

• Stochastic Latent World Modeling: We propose a stochastic latent transition framework that generates multi-step rollouts in representation space using Gaussian parameterization and recurrent refinement, enabling predictive structure without requiring temporal supervision.

• Latent Rollout Planning Mechanism: We introduce an adaptive rollout aggregation strategy that learns to weight latent states, producing a trajectory-aware representation for classification.

• Structured Multi-Objective Learning: We design a unified training objective combining masked reconstruction, spectral and spatial latent dynamics, stochastic rollout regularization, and physically motivated constraints, ensuring both discriminative power and structural consistency.

The remainder of this paper is organized as follows. Section II introduces the proposed Stochastic Latent World Modeling (LWM) framework and its learning objectives. Section III details the experimental setup. Sections IV-V analyze latent rollout reasoning and robustness, while Section VI examines model interpretability. Section VII presents ablation studies, and Section VIII compares LWM with state-of-the-art methods. Section IX discusses limitations and failure cases. Finally, Section X concludes the paper and outlines future directions.

Mathematical Modeling

The proposed LWM model formulates HSIC as a trajectory-aware spatial-spectral representation learning problem. Each input patch is decomposed into spectral and spatial token sequences and encoded through dedicated Transformer branches, whose outputs are fused into a joint representation. From this representation, a compact environment state is constructed and propagated through a stochastic LWM module to generate multi-step rollouts in representation space. Importantly, the term “world modeling” does not imply temporal simulation: the latent state represents the observed spatial-spectral context, while the learned stochastic transition operator maps this state to plausible future latent configurations under the structural dynamics captured by the model. Accordingly, a “state transition” denotes a learned transformation in representation space, and a “rollout” denotes repeated application of this transition to generate a sequence of latent hypotheses. These latent trajectories do not correspond to explicit temporal observations or imply that the hyperspectral scene evolves over time.

For static HSIC, the purpose of such latent transitions is to provide structurally grounded alternatives to unavailable temporal supervision. Specifically, spectral adjacency, spatial neighborhood relations, and their fused context define the relational structure of the latent transition space. The resulting rollout therefore explores plausible representations consistent with the observed patch structure, rather than simulating physical or temporal evolution. A rollout planner aggregates the latent trajectory into a planned state, which is combined with the current state to produce trajectory-aware class predictions. In this way, the trajectory serves both as a predictive regularizer and as a decision mechanism, allowing the classifier to aggregate multiple structurally plausible latent alternatives rather than relying exclusively on a single static embedding. The overall procedure is summarized in Algorithm 1 and Figure 1.

Problem Formulation and Notation

HSIC assigns a semantic label to each spatial location using spectral signatures and local spatial context. Let the HSI be $\mathbf{X} \in\mathbb{R}^{H_0 \times W_0 \times C}$, where $H_0$, $W_0$, and $C$ denote spatial dimensions and spectral bands. For a labeled pixel at spatial coordinates $\mathbf{q}_i=(r_i,c_i)$, where $r_i$ and $c_i$ denote its row and column indices, respectively, a local odd patch of size $P \times P$ is extracted using mirror padding, yielding $\mathbf{x}_i \in\mathbb{R}^{C \times P \times P}$. For a mini-batch of $B$ samples, the input is $\mathbf{x} \in\mathbb{R}^{B \times C \times P \times P}$ with labels $\mathbf{y} \in\{0,\ldots,K-1\}^{B}$, where $K$ is the number of classes.

Spectral and Spatial Tokenization

The first step transforms each hyperspectral patch into complementary spectral and spatial token representations [14].

Spectral Tokenization. The spectral tokenizer constructs tokens along the wavelength axis, where each spectral band is treated as an individual token retaining its local spatial context. For the $i$-th sample, the $c$-th spectral band $\mathbf{x}_{i,c,:,:} \in\mathbb{R}^{P \times P}$ is flattened into a vector $\mathbf{u}_{i,c} \in\mathbb{R}^{P^2}$ and projected into a $D$-dimensional embedding space. A learnable spectral positional encoding is added to preserve the ordering of wavelength indices. For a mini-batch, this process is expressed as:

where $\mathbf{x}_{\mathrm{spec}} \in\mathbb{R}^{B \times C \times P^{2}}$ is the flattened spectral representation, $\mathbf{W}_{\mathrm{spec}} \in\mathbb{R}^{P^{2} \times D}$ is a learnable projection matrix, and $\mathbf{P}_{\lambda} \in\mathbb{R}^{1 \times C \times D}$ denotes the learnable spectral positional embedding. The batch-expanded positional embedding is defined as $\widetilde{\mathbf{P}}_{\lambda} = \mathbf{1}_{B} \otimes\mathbf{P}_{\lambda} \in\mathbb{R}^{B \times C \times D}$, where $\mathbf{1}_{B}$ is a $B$-dimensional vector of ones and $\otimes$ denotes replication along the batch dimension.

Spatial Tokenization. The spatial tokenizer constructs tokens along the spatial domain, where each token corresponds to a spatial location and carries its full spectral signature. The input tensor $\mathbf{x} \in\mathbb{R}^{B \times C \times P \times P}$ is rearranged by flattening the two spatial dimensions into a single spatial-token dimension, such that $\mathbf{x}_{\mathrm{spat}} \in\mathbb{R}^{B \times N \times C}$ with $N = P^{2}$. Thus, for the $i$-th sample and the $p$-th spatial location, the corresponding spectral vector is $\mathbf{x}_{\mathrm{spat},i,p,:} \in\mathbb{R}^{C}$, which is projected into the same $D$-dimensional embedding space. A learnable spatial positional encoding is added to retain spatial layout information. The resulting spatial token sequence is given by

where $\mathbf{x}_{\mathrm{spat}} \in\mathbb{R}^{B \times N \times C}$ is the spatially rearranged tensor, $\mathbf{W}_{\mathrm{spat}} \in\mathbb{R}^{C \times D}$ is a learnable projection matrix, and $\mathbf{P}_{xy} \in\mathbb{R}^{1 \times N \times D}$ represents the spatial positional embeddings. The batch-expanded positional embedding is defined as $\widetilde{\mathbf{P}}_{xy} = \mathbf{1}_{B} \otimes\mathbf{P}_{xy} \in\mathbb{R}^{B \times N \times D}$.

Transformer Encoding

Let $\operatorname{SpecEnc}_{\Theta_{\lambda}}(\cdot)$ denote the spectral transformer encoder composed of $L_{\lambda}$ transformer layers, and let $\operatorname{SpatEnc}_{\Theta_{s}}(\cdot)$ denote the spatial transformer encoder composed of $L_{s}$ transformer layers. The corresponding encoded representations are defined as

where $C$ and $N = P^{2}$ denote the sequence lengths of the spectral and spatial branches, respectively. The parameters $\Theta_{\lambda}$ and $\Theta_{s}$ are learned independently. For a token sequence $\mathbf{T}^{\ell-1} \in\mathbb{R}^{B \times M \times D}$ at layer $\ell- 1$, a transformer block updates the representation through pre-normalized self-attention and feed-forward sublayers:

where $\operatorname{MHSA}(\cdot)$ denotes multi-head self-attention [17], $\operatorname{MLP}(\cdot)$ denotes a feed-forward network, $\operatorname{LN}(\cdot)$ denotes layer normalization, and $\operatorname{DropPath}(\cdot)$ introduces stochastic depth. The sequence length $M$ is set to $C$ for the spectral branch and $N$ for the spatial branch. Applying this transformer block for $L_{\lambda}$ and $L_{s}$ layers yields the encoded representations $\mathbf{S}$ and $\mathbf{A}$, respectively.

Fusion. To enable joint modeling, the spectral and spatial representations are fused using a transformer encoder (using $L_{f}$ transformer layers). A learnable classification token $\mathbf{c}_{0} \in\mathbb{R}^{1 \times1 \times D}$ and type embeddings $\mathbf{e}_{\mathrm{cls}}, \mathbf{e}_{\lambda}, \mathbf{e}_{s} \in\mathbb{R}^{1 \times1 \times D}$ are introduced to distinguish token types. Since these embeddings are shared across the mini-batch, their batch-expanded forms are defined as $\widetilde{\mathbf{c}}_{0} = \mathbf{1}_{B} \otimes\mathbf{c}_{0} \in\mathbb{R}^{B \times1 \times D}$, $\widetilde{\mathbf{e}}_{\mathrm{cls}} = \mathbf{1}_{B} \otimes\mathbf{e}_{\mathrm{cls}} \in\mathbb{R}^{B \times1 \times D}$, $\widetilde{\mathbf{e}}_{\lambda} = \mathbf{1}_{B} \otimes\mathbf{e}_{\lambda} \in\mathbb{R}^{B \times1 \times D}$, and $\widetilde{\mathbf{e}}_{s} = \mathbf{1}_{B} \otimes\mathbf{e}_{s} \in\mathbb{R}^{B \times1 \times D}$. Let $\operatorname{Broadcast}_{C}(\cdot)$ and $\operatorname{Broadcast}_{N}(\cdot)$ denote replication along the spectral and spatial sequence dimensions, respectively, and let $\operatorname{Concat}_{\mathrm{seq}}(\cdot)$ denote concatenation along the sequence dimension. The fused input sequence is defined as

The fused sequence is processed by the fusion transformer encoder as $\mathbf{F} = \operatorname{Fuse}_{\Theta}(\mathbf{F}^{0}) \in\mathbb{R}^{B \times(1+C+N) \times D}$. The output is decomposed as $\mathbf{h}_{\mathrm{cls}} = \mathbf{F}_{:,0,:} \in\mathbb{R}^{B \times D}$, $\mathbf{Z}_{\lambda} = \mathbf{F}_{:,1:1+C,:} \in\mathbb{R}^{B \times C \times D}$, and $\mathbf{Z}_{s} = \mathbf{F}_{:,1+C:,:} \in\mathbb{R}^{B \times N \times D}$. The classification token $\mathbf{h}_{\mathrm{cls}}$ provides a global representation for the base prediction. The base classification head $g_{\mathrm{base}} : \mathbb{R}^{D} \to\mathbb{R}^{K}$ consists of layer normalization, a linear layer, GELU activation, dropout, and a final linear layer mapping to $K$ classes. The base logits are given by

and serve as an initial prediction that is refined through latent rollout modeling in subsequent stages.

Environment State Construction

Following spatial-spectral encoding and fusion, a compact latent representation is constructed to capture the local environment of the target pixel. While the classification token $h_{\mathrm{cls}}$ provides a global summary, it does not explicitly separate local spatial structure and spectral composition. To address this, the environment state is formed by aggregating complementary components derived from the fused tokens. Let $n_0 = \frac{P^2-1}{2}$ denote the flattened index of the center location in the spatial token sequence, assuming $P$ is odd, and let $\mathcal{N}_8(n_0)$ denote the set of indices corresponding to its eight first-ring neighbors. The center state is defined as $z_c = Z_{s,:,n_0} \in\mathbb{R}^{B \times N \times D}$, which encodes the latent representation of the target pixel. The neighbor state is computed as:

capturing local spatial context. The global spectral state is given by:

and the scene-level representation is $z_g = h_{\mathrm{cls}} \in\mathbb{R}^{B \times D}$. These components are concatenated along the column dimension using the operator $\operatorname{Concat}_{\mathrm{col}}(\cdot)$, yielding $z_{\mathrm{cat}} = \operatorname{Concat}_{\mathrm{col}}(z_c, z_n, z_\lambda, z_g) \in\mathbb{R}^{B \times4D}$. The concatenated representation is then mapped into a unified latent representation:

where $E_{\mathrm{env}} : \mathbb{R}^{B \times4D} \to\mathbb{R}^{B \times D}$ is implemented using layer normalization, linear projection, Gaussian error linear unit (GELU) activation, dropout, and a final normalization layer. The resulting $z_0$ serves as the initial latent environment state, integrating local spatial, spectral, and global contextual information for subsequent latent rollout modeling.

Stochastic Latent World Modeling

Given the environment state $z_0$, the model simulates its evolution in latent space through a stochastic rollout mechanism [42, 43]. This process captures variations of the local spatial-spectral representation by modeling uncertainty in latent transitions [50]. The stochastic LWM module defines a multi-step rollout starting from $z_0$ [44]. For $t = \{0, \ldots, T-1\}$, a transition network $G_\theta$ predicts the parameters of a Gaussian latent transition:

where $G_\theta: \mathbb{R}^{B \times D} \to\mathbb{R}^{B \times2D}$ is implemented as a feed-forward network consisting of layer normalization, a linear layer mapping $D$ to $2D$, GELU activation, dropout, and a final linear layer mapping $2D$ to $2D$. The output is split along the feature dimension into the predicted mean $\mu_{t+1} \in\mathbb{R}^{B \times D}$ and log-variance $\log\sigma_{t+1}^{2} \in\mathbb{R}^{B \times D}$. For numerical stability, the predicted log-variance is clipped element-wise to a fixed interval $[a,b]$, where $a < b$ are predefined lower and upper bounds, respectively:

where, for each element $v$ of $\log\sigma_{t+1}^{2}$ it is $\operatorname{clip}(v,a,b) = \min\{\max\{v,a\},b\}$. Thus, all predicted log-variance values are constrained to the interval $[a,b]$. During training, the next latent state is sampled using the reparameterization trick [39, 51]:

where $\epsilon_{t+1} \sim\mathcal{N}_{B \times D}(0_{B \times D}, I_B \otimes I_D)$, so that $\epsilon_{t+1} \in\mathbb{R}^{B \times D}$ and $\odot$ denotes element-wise Hadamard multiplication. During inference, stochastic sampling is disabled and $\widetilde{z}_{t+1} = \mu_{t+1}$. The state is then refined via the recurrent update

where $R_\theta$ is implemented as a gated recurrent unit (GRU) cell [45]. This yields a latent trajectory. Since each latent state, predicted mean, and log-variance is represented as a matrix in $\mathbb{R}^{B \times D}$, we define $\operatorname{Concat}_3(\cdot)$ as concatenation along a new third dimension. For $T$ matrices $\{X_1,\ldots,X_T\} \in\mathbb{R}^{B \times D}$ as $\operatorname{Concat}_3(X_1,\ldots,X_T) \in\mathbb{R}^{B \times D \times T}$, and $\operatorname{reshape}_{B \times T \times D}(\cdot)$ denotes the corresponding reshaping of this into the batch-trajectory-feature ordering. The latent trajectory is therefore defined as $Z_{\mathrm{roll}} = \operatorname{reshape}_{B \times(T+1) \times D}(\operatorname{Concat}_3(z_0,z_1,\ldots,z_T))$, where $Z_{\mathrm{roll}} \in\mathbb{R}^{B \times(T+1) \times D}$.

Similarly, the predicted means and log-variances are collected as $\boldsymbol{\mu} = \operatorname{reshape}_{B \times T \times D}(\operatorname{Concat}_3(\mu_1,\ldots,\mu_T)) \in\mathbb{R}^{B \times T \times D}$, and $\mathbf{V} = \operatorname{reshape}_{B \times T \times D}(\operatorname{Concat}_3(\log\sigma_1^2,\ldots,\log\sigma_T^2)) \in\mathbb{R}^{B \times T \times D}$.

To construct an auxiliary latent target, let $Q_\theta: \mathbb{R}^{B \times D} \to\mathbb{R}^{B \times D}$ denote a learnable latent projection network that maps the current latent state to a target representation. The stop-gradient operator $\operatorname{sg}(\cdot)$ returns its input while preventing gradients from propagating through that input during back-propagation, i.e., $\operatorname{sg}(u) = u$, $\frac{\partial\operatorname{sg}(u)}{\partial u} = 0$. The auxiliary latent target is then defined as:

and the corresponding sequence is constructed using the same third-dimension concatenation and reshaping operation $\widehat{Z} = \operatorname{reshape}_{B \times T \times D}(\operatorname{Concat}_3(\widehat{z}_1,\ldots,\widehat{z}_T)) \in\mathbb{R}^{B \times T \times D}$. This target acts as an internal projection used to regularize the stochastic transitions. The rollout operates entirely in representation space and does not correspond to explicit temporal observations. The resulting latent trajectory provides a set of states that are subsequently used for trajectory-aware aggregation and classification.

Trajectory Planning and Classification

Given the latent trajectory $\mathbf{z}_{\{t=0\}}^{T}$, the model performs trajectory-aware aggregation to obtain a decision-relevant representation. For each rollout state $\mathbf{z}_t \in\mathbb{R}^{B \times D}$, a learnable scoring network $u_\theta: \mathbb{R}^{D} \to\mathbb{R}$ assigns a scalar score independently to each sample. Specifically, $u_\theta$ is implemented as a two-layer MLP with layer normalization, a linear projection from $D$ to $D/2$, GELU activation, and a linear projection to a scalar:

where $\mathbf{W}_u \in\mathbb{R}^{D/2 \times D}$, $\mathbf{b}_u \in\mathbb{R}^{D/2}$, $\mathbf{w}_u \in\mathbb{R}^{D/2}$, and $b_u \in\mathbb{R}$ are learnable parameters. Thus, for sample $i$ and rollout step $t$, the corresponding score is $\alpha_{i,t} = u_\theta(\mathbf{z}_{i,t}) \in\mathbb{R}$, $i = \{1,\ldots,B\}$, $t = \{0,\ldots,T\}$, where $\mathbf{z}_{i,t} \in\mathbb{R}^{D}$ denotes the $i$-th row of $\mathbf{z}_t$. The scores are normalized independently for each sample along the trajectory dimension using the softmax operator:

where $\sum_{t=0}^{T} w_{i,t} = 1$, $w_{i,t} \ge0$ and $\boldsymbol{\alpha}_i = [\alpha_{i,0},\ldots,\alpha_{i,T}] \in\mathbb{R}^{T+1}$. Equivalently, $\operatorname{softmax}(\boldsymbol{\alpha}_i)$ converts the scalar scores associated with the $T+1$ rollout states of sample $i$ into normalized trajectory weights. The planned representation is then obtained by a weighted aggregation of the latent states $\mathbf{z}_{\mathrm{plan},i} = \sum_{t=0}^{T} w_{i,t}\mathbf{z}_{i,t} \in\mathbb{R}^{D}$, for the complete mini-batch,

where $\mathbf{w}_t = [w_{1,t},\ldots,w_{B,t}]^{\mathsf{T}} \in\mathbb{R}^{B}$ and the weights are broadcast across the feature dimension. The trajectory representation is constructed by concatenating the initial latent state and the planned state along the feature dimension:

The trajectory-based logits are produced by a trajectory classifier $g_{\mathrm{traj}} : \mathbb{R}^{2D} \to\mathbb{R}^{K}$ which consists of layer normalization, a linear projection from $2D$ to $D$, GELU activation, dropout, and a final linear projection from $D$ to $K$ as $g_{\mathrm{traj}}(\mathbf{r}) = \mathbf{W}_{\mathrm{traj},2}\operatorname{GELU}(\mathbf{W}_{\mathrm{traj},1}\operatorname{LN}(\mathbf{r}) + \mathbf{b}_{\mathrm{traj},1})+\mathbf{b}_{\mathrm{traj},2}$, where $\mathbf{W}_{\mathrm{traj},1} \in\mathbb{R}^{D \times2D}$, $\mathbf{b}_{\mathrm{traj},1} \in\mathbb{R}^{D}$, $\mathbf{W}_{\mathrm{traj},2} \in\mathbb{R}^{K \times D}$, and $\mathbf{b}_{\mathrm{traj},2} \in\mathbb{R}^{K}$ are learnable parameters, and dropout is applied between the GELU activation and the linear projection. Thus,

The final prediction combines the base and trajectory logits using a sample-dependent branch gate. The gate is computed from the concatenation of the fused CLS representation and the planned latent state as:

where $\mathbf{h}_{\mathrm{cls},i} \in\mathbb{R}^{D}$ and $\mathbf{z}_{\mathrm{plan},i} \in\mathbb{R}^{D}$ denote the $i$-th samples of the CLS and planned representations, respectively. $\mathbf{W}_g \in\mathbb{R}^{D \times2D}$, $\mathbf{b}_g \in\mathbb{R}^{D}$, $\mathbf{w}_g \in\mathbb{R}^{D}$, and $b_g \in\mathbb{R}$ are learnable parameters, and $\sigma(\cdot)$ denotes the sigmoid function. The stop-gradient operation $\operatorname{sg}(\cdot)$ prevents gradients from the branch-gating network from propagating through the planned state. Collecting the sample-wise gates gives $\mathbf{g} = [g_1,\ldots,g_B]^{\mathsf{T}} \in(0,1)^{B \times1}$. To apply the scalar gate to all $K$ class logits of each sample, define $\mathbf{G} = \mathbf{g}\mathbf{1}_K^{\mathsf{T}} \in(0,1)^{B \times K}$, where $\mathbf{1}_K \in\mathbb{R}^{K}$ is a vector of ones. The final logits are then obtained by element-wise interpolation:

where $\odot$ denotes element-wise multiplication and $\mathbf{1}_{B \times K}$ is a $B \times K$ matrix of ones. Thus, each sample uses a single learned gate to control the contribution of the base and trajectory branches across all classes. For the predicted label, let $\boldsymbol{\ell}_i \in\mathbb{R}^{K}$ denote the $i$-th row of the final logit matrix $\boldsymbol{\ell} \in\mathbb{R}^{B \times K}$. The predicted class is given by

This trajectory planning mechanism enables the model to aggregate information across the latent rollout while retaining the initial environment state for trajectory-based classification. By learning sample-dependent importance weights over the rollout states, the model incorporates structured information from the latent trajectory into the final prediction.

Reconstruction and Dynamics Constraints

The model incorporates reconstruction and latent dynamics constraints to regularize the spectral and spatial token representations. These constraints operate on the fused spectral and spatial representations $\mathbf{Z}_{\lambda}$ and $\mathbf{Z}_s$ obtained from the fusion transformer. In particular, the reconstruction heads recover the corresponding input-domain representations, while the dynamics modules impose predictive structure along the spectral sequence and across the local spatial neighborhood.

Spectral Reconstruction. The spectral reconstruction head $R_{\lambda}$ maps each $D$-dimensional spectral token back to the flattened spatial patch associated with the corresponding spectral band. Given $\mathbf{Z}_{\lambda} \in\mathbb{R}^{B \times C \times D}$, the same reconstruction mapping is applied independently to all $C$ spectral tokens. Specifically, $R_{\lambda} : \mathbb{R}^{D} \to\mathbb{R}^{N}$ is implemented using layer normalization, a linear projection from $D$ to $D$, GELU activation, and a final linear projection from $D$ to $N$:

where $\mathbf{W}^{\mathrm{rec}}_{\lambda,1} \in\mathbb{R}^{D \times D}$, $\mathbf{b}^{\mathrm{rec}}_{\lambda,1} \in\mathbb{R}^{D}$, $\mathbf{W}^{\mathrm{rec}}_{\lambda,2} \in\mathbb{R}^{N \times D}$, and $\mathbf{b}^{\mathrm{rec}}_{\lambda,2} \in\mathbb{R}^{N}$ are learnable parameters. Applying $R_{\lambda}$ to every spectral token gives $\hat{\mathbf{x}}_{\lambda} = R_{\lambda}(\mathbf{Z}_{\lambda}) \in\mathbb{R}^{B \times C \times N}$. The corresponding target is obtained by flattening the two spatial dimensions of the input patch $\mathbf{x}_{\lambda} = \operatorname{reshape}(\mathbf{x}) \in\mathbb{R}^{B \times C \times N}$. Thus, the $c$-th output token of $R_{\lambda}$ reconstructs the $P \times P$ spatial values belonging to the $c$-th spectral band.

Spatial Reconstruction. The spatial reconstruction head $R_s$ maps each spatial token back to its corresponding spectral vector. Given $\mathbf{Z}_s \in\mathbb{R}^{B \times N \times D}$, the same reconstruction mapping is applied independently to all $N$ spatial tokens. The mapping $R_s : \mathbb{R}^{D} \to\mathbb{R}^{C}$ is implemented using layer normalization, a linear projection, GELU activation, and a linear projection from $D$ to $C$:

where $\mathbf{W}^{\mathrm{rec}}_{s,1} \in\mathbb{R}^{D \times D}$, $\mathbf{b}^{\mathrm{rec}}_{s,1} \in\mathbb{R}^{D}$, $\mathbf{W}^{\mathrm{rec}}_{s,2} \in\mathbb{R}^{C \times D}$, and $\mathbf{b}^{\mathrm{rec}}_{s,2} \in\mathbb{R}^{C}$ are learnable parameters. Applying $R_s$ to all spatial tokens yields $\hat{\mathbf{x}}_s = R_s(\mathbf{Z}_s) \in\mathbb{R}^{B \times N \times C}$. The corresponding target is obtained by first permuting the input tensor so that the spatial dimensions precede the spectral dimension and then flattening the spatial dimensions $\mathbf{x}_s = \operatorname{reshape}(\operatorname{permute}(\mathbf{x})) \in\mathbb{R}^{B \times N \times C}$. Hence, each spatial token is trained to reconstruct the complete $C$-dimensional spectral signature of its corresponding spatial location. When masking is applied, the reconstruction losses are evaluated only on the masked spectral bands or spatial locations; otherwise, all entries are used. Let $\rho_\beta(\cdot,\cdot)$ denote the element-wise Smooth-$\ell_1$ penalty with transition parameter $\beta> 0$:

The spectral and spatial reconstruction losses are then defined as $\mathcal{L}_{\mathrm{spec\text{-}rec}} = \operatorname{mean}\left[\rho_{\beta_{\mathrm{rec}}}(\hat{\mathbf{x}}_\lambda,\mathbf{x}_\lambda)\right]$ and $\mathcal{L}_{\mathrm{spat\text{-}rec}} = \operatorname{mean}\left[\rho_{\beta_{\mathrm{rec}}}(\hat{\mathbf{x}}_s,\mathbf{x}_s)\right]$, where the mean is taken over the selected masked entries when masking is enabled and over all entries otherwise.

Spectral Dynamics. The spectral dynamics module models the transition from one spectral token to the next. Given $\mathbf{Z}_\lambda= [\mathbf{z}_{\lambda,1},\ldots,\mathbf{z}_{\lambda,C}]$ with $\mathbf{z}_{\lambda,c} \in\mathbb{R}^{B \times D}$, the initial recurrent state is the first spectral token $\mathbf{h}_{\lambda,1} = \mathbf{Z}_{\lambda,:,1,:} \in\mathbb{R}^{B \times D}$. For each subsequent spectral position $c = \{2,\ldots,C\}$, the recurrent state is updated using a GRU cell:

The updated state is mapped to the predicted representation of the next spectral token using a prediction head consisting of layer normalization, a linear projection, GELU activation, dropout, and a linear projection from $D$ to $D$:

where $H_\lambda(\mathbf{h}) = \mathbf{W}^{\mathrm{dyn}}_{\lambda,2}\operatorname{Dropout}\left[\operatorname{GELU}\left(\mathbf{W}^{\mathrm{dyn}}_{\lambda,1}\operatorname{LN}(\mathbf{h}) + \mathbf{b}^{\mathrm{dyn}}_{\lambda,1}\right)\right] + \mathbf{b}^{\mathrm{dyn}}_{\lambda,2}$. The resulting predictions are collected as $\hat{\mathbf{Z}}_\lambda^{+} = \operatorname{stack}(\hat{\mathbf{z}}_{\lambda,2},\ldots,\hat{\mathbf{z}}_{\lambda,C}) \in\mathbb{R}^{B \times(C-1) \times D}$. The corresponding target is the sequence of subsequent spectral tokens with gradients detached:

Thus, $\hat{\mathbf{Z}}_\lambda^{+}$ is not a slice of the latent rollout tensor $\hat{\mathbf{Z}}$; it is the sequence of next-band predictions generated directly from the spectral tokens by the spectral dynamics module. The spectral dynamics loss is defined as the mean-squared error between the predicted next-band representations and their detached targets $\mathcal{L}_{\mathrm{spec\text{-}dyn}} = \operatorname{mean}\left[(\hat{\mathbf{Z}}_\lambda^{+} - \mathbf{Z}_\lambda^{+})^2\right]$.

Spatial Dynamics. The spatial dynamics module models the relationship between the center spatial token and its eight first-ring neighbors. Let $\mathcal{N}_8(n_0)$ denote the set of the eight spatial indices obtained from the offsets $(\Delta_r,\Delta_c) \in\{-1,0,1\}^2 \setminus(0,0)$, around the center location. The corresponding neighboring spatial tokens are $\mathbf{Z}_{s,\mathcal{N}} = \mathbf{Z}_{s,:, \mathcal{N}_8(n_0),:} \in\mathbb{R}^{B \times8 \times D}$, while the center token is $\mathbf{z}_{s,c} = \mathbf{Z}_{s,:,n_0,:} \in\mathbb{R}^{B \times D}$. The center token is passed through a prediction head consisting of layer normalization, a linear projection, GELU activation, dropout, and a linear projection from $D$ to $D$:

where $H_s(\mathbf{z}) = \mathbf{W}^{\mathrm{dyn}}_{s,2}\operatorname{Dropout}\left[\operatorname{GELU}\left(\mathbf{W}^{\mathrm{dyn}}_{s,1}\operatorname{LN}(\mathbf{z}) + \mathbf{b}^{\mathrm{dyn}}_{s,1}\right)\right] + \mathbf{b}^{\mathrm{dyn}}_{s,2}$. The same predicted representation is compared with each of the eight first-ring neighbors $\hat{\mathbf{Z}}_{s,\mathcal{N}} = \operatorname{Broadcast}_8(\hat{\mathbf{z}}_{s,n}) \in\mathbb{R}^{B \times8 \times D}$, where $\operatorname{Broadcast}_8(\cdot)$ replicates the predicted center representation across each neighbor positions. The corresponding target is detached from the encoder $\mathbf{Z}^{\mathrm{stop}}_{s,\mathcal{N}} = \operatorname{sg}(\mathbf{Z}_{s,:, \mathcal{N}_8(n_0),:}) \in\mathbb{R}^{B \times8 \times D}$. The spatial dynamics loss is therefore $\mathcal{L}_{\mathrm{spat\text{-}dyn}} = \operatorname{mean}\left[(\hat{\mathbf{Z}}_{s,\mathcal{N}}-\mathbf{Z}^{\mathrm{stop}}_{s,\mathcal{N}})^2\right]$.

Composite Training Objective. The complete training objective combines the final classification loss, base classification loss, reconstruction constraints, spectral and spatial dynamics losses, stochastic latent rollout regularization, trajectory supervision, trajectory consistency, and the additional spectral and spatial structural constraints. To avoid fixed weighting constants, all relative contributions are controlled by nonnegative learnable training hyperparameters:

where all coefficients $\lambda_{\mathrm{final}}, \lambda_{\mathrm{base}}, \lambda_{\mathrm{sr}}, \lambda_{\mathrm{pr}}, \lambda_{\mathrm{sd}}, \lambda_{\mathrm{pd}}, \lambda_{\mathrm{env}}, \lambda_{\mathrm{KL}}, \lambda_{\mathrm{traj}}, \lambda_{\mathrm{tc}}, \lambda_{\mathrm{SAM}}, \lambda_{\mathrm{sm}}, \lambda_{\mathrm{sc}} \ge0$ are hyperparameters. The final and base classification losses use label-smoothed cross-entropy. Let $\ell_i \in\mathbb{R}^{K}$ denote the final logits for sample $i$, $\ell_{\mathrm{base},i} \in\mathbb{R}^{K}$ the corresponding base logits, and $y_i \in\{0,\ldots,K-1\}$ the ground-truth (GT) class. With label-smoothing parameter $\epsilon_{\mathrm{ls}} \in[0,1)$, define $q_{i,k} = (1-\epsilon_{\mathrm{ls}})[k=y_i] + \frac{\epsilon_{\mathrm{ls}}}{K}$ and $\mathcal{L}_{\mathrm{final}} = -\frac{1}{B}\sum_{i=1}^{B}\sum_{k=0}^{K-1}q_{i,k}\log\operatorname{softmax}(\ell_i)_k$ and $\mathcal{L}_{\mathrm{base}} = -\frac{1}{B}\sum_{i=1}^{B}\sum_{k=0}^{K-1}q_{i,k}\log\operatorname{softmax}(\ell_{\mathrm{base},i})_k$.

The trajectory-classification loss $\mathcal{L}_{\mathrm{traj\text{-}cls}}$ is defined analogously using the trajectory logits $\ell_{\mathrm{traj},i}$ as $\mathcal{L}_{\mathrm{traj\text{-}cls}} = -\frac{1}{B}\sum_{i=1}^{B}\sum_{k=0}^{K-1}q_{i,k}\log\operatorname{softmax}(\ell_{\mathrm{traj},i})_k$.

The trajectory-consistency loss encourages the trajectory branch to remain consistent with the detached base prediction. Let $\mathbf{p}_{\mathrm{base}} = \operatorname{softmax}(\operatorname{sg}(\ell_{\mathrm{base}}))$ and $\log\mathbf{p}_{\mathrm{traj}} = \log\operatorname{softmax}(\ell_{\mathrm{traj}})$. Then $\mathcal{L}_{\mathrm{traj\text{-}cons}} = \frac{1}{B}\sum_{i=1}^{B}\sum_{k=0}^{K-1}p_{\mathrm{base},i,k}\log\frac{p_{\mathrm{base},i,k}}{p_{\mathrm{traj},i,k}}$. For the stochastic latent rollout, let $\widehat{\mathbf{Z}} \in\mathbb{R}^{B \times T \times D}$ denote the detached rollout targets, let $\boldsymbol{\mu} \in\mathbb{R}^{B \times T \times D}$ denote the predicted Gaussian means, and let $\mathbf{V} \in\mathbb{R}^{B \times T \times D}$ denote the predicted log-variances. The Gaussian negative log-likelihood is

The rollout smoothness term penalizes abrupt changes between consecutive latent states $\mathcal{L}_{\mathrm{smooth-roll}} = \operatorname{mean}\left[\left|\mathbf{Z}_{\mathrm{roll}}[:,1:,:]-\mathbf{Z}_{\mathrm{roll}}[:,:-1,:]\right|\right]$, where $\mathbf{Z}_{\mathrm{roll}} \in\mathbb{R}^{B \times(T+1) \times D}$ denotes the latent rollout trajectory. The complete environment rollout loss is $\mathcal{L}_{\mathrm{env}} = \mathcal{L}_{\mathrm{NLL}} + \lambda_{\mathrm{smooth-roll}}\mathcal{L}_{\mathrm{smooth-roll}}$, where $\lambda_{\mathrm{smooth-roll}} \ge0$ is a hyperparameter rather than a fixed constant. The latent Gaussian regularization term is the KL divergence between the predicted Gaussian transition distribution and a standard normal prior $\mathcal{L}_{\mathrm{KL}} = -\frac{1}{2}\operatorname{mean}\left[1+\mathbf{V}-\boldsymbol{\mu}^2-\exp(\mathbf{V})\right]$.

Finally, the spectral-angle, spectral-smoothness, and spatial-consistency terms are defined from the corresponding model outputs. Let $\widehat{\mathbf{x}}_{\lambda}$ denote the spectral reconstruction. The spectral-angle loss is computed on the center spatial position of each reconstructed patch:

where $\widehat{\mathbf{v}}_i,\mathbf{v}_i \in\mathbb{R}^{C}$ are respectively the reconstructed and true center-pixel spectral vectors and $\epsilon_{\mathrm{SAM}} > 0$ is a numerical-stability parameter. The spectral smoothness loss penalizes changes between adjacent reconstructed spectral bands $\mathcal{L}_{\mathrm{smooth}} = \operatorname{mean}\left[\left|\widehat{\mathbf{x}}_{\lambda}[:,2,:]-\widehat{\mathbf{x}}_{\lambda}[:,1,:]\right|\right]$.

For a batch containing spatially neighboring labeled pixels, let $\mathbf{p}_i=\operatorname{softmax}(\boldsymbol{\ell}_i)$ and let $\mathcal{P}$ denote the set of sample pairs whose spatial coordinates are within the prescribed local neighborhood. The spatial-consistency loss is $\mathcal{L}_{\mathrm{spat-cons}} = \frac{1}{|\mathcal{P}|}\sum_{(i,j)\in\mathcal{P}}\lVert\mathbf{p}_i-\mathbf{p}_j\rVert^2$. When no valid neighboring pair is present in the mini-batch, $\mathcal{L}_{\mathrm{spat-cons}}$ is set to zero.

Experimental Setup

The proposed LWM framework evaluates trajectory-aware spatial-spectral representation learning for HSIC in a within-scene, label-scarce setting, where limited labeled samples from a fixed hyperspectral scene are used to learn structured latent representations for classifying held-out samples from the same scene.

Datasets. We evaluate the proposed LWM framework on five publicly available HSI benchmarks: WHU-HongHu (HH) [52], Tangdaowan (TD) [53], DFC18 (training set only) [54], a subset of HyRANK e.g., Dioni, and Loukia [55]. These datasets encompass a broad range of scene types, spatial resolutions, spectral characteristics, and classification complexities, providing a comprehensive evaluation of the proposed framework across agricultural, urban, coastal, and mixed land-cover environments.

Implementation Details. The proposed LWM model is optimized using AdamW with an initial learning rate of $2 \times10^{-4}$ and weight decay of $7 \times10^{-4}$. A cosine annealing learning-rate scheduler is employed over 100 training epochs, while early stopping with a patience of 15 epochs is used to prevent overfitting. Gradient clipping with a maximum norm of 1.0 is applied for training stability, and mixed-precision training is utilized to improve computational efficiency. For all datasets, hyperspectral cubes are normalized using per-band standardization. Local spatial neighborhoods are extracted using mirror padding with a fixed patch size of $25 \times25$. Unless otherwise specified, the labeled samples are divided into training, validation, and testing subsets using stratified sampling, with ratios of 2%, 2%, and 96% for the Dioni, HH, and TD datasets, and 10%, 10%, and 80% for the DFC18 and Loukia datasets, respectively. Each selected pixel serves as the center of its corresponding spatial-spectral patch; thus, the sampling ratios refer to labeled center pixels. Unless otherwise stated, a batch size of 64 is used, and label smoothing is set to 0.1.

The transformer backbone employs an embedding dimension of $D=128$, with spectral, spatial, and fusion transformer depths set to 2. Each transformer layer uses 8 attention heads, an MLP expansion ratio of 4.0, dropout of 0.1, and stochastic depth with a maximum drop-path rate of 0.05. During training, random spectral-band masking and spatial-pixel masking are applied with masking ratios of 25% and 15%, respectively. For LWM module, the rollout horizon is fixed to $T=3$. The stochastic transition model predicts Gaussian latent distributions with log-variance constrained to the interval $[a,b]$, where the lower and upper bounds are fixed as $a=-6$ and $b=2$. Rollout states are refined through a GRU-based latent transition module and aggregated using the proposed rollout planner. The latent state dimension is set equal to the transformer embedding dimension ($D=128$). The overall objective combines classification, reconstruction, latent dynamics, rollout regularization, trajectory supervision, and physically motivated constraints. Unless otherwise specified, the corresponding weighting coefficients are set to $\lambda_{\mathrm{spec-rec}}=0.25$, $\lambda_{\mathrm{spat-rec}}=0.25$, $\lambda_{\mathrm{spec-dyn}}=0.20$, $\lambda_{\mathrm{spat-dyn}}=0.20$, $\lambda_{\mathrm{roll}}=0.25$, $\lambda_{\mathrm{KL}}=0.001$, $\lambda_{\mathrm{traj-cons}}=0.10$, $\lambda_{\mathrm{traj-cls}}=0.25$, $\lambda_{\mathrm{SAM}}=0.08$, $\lambda_{\mathrm{smooth}}=0.05$, and $\lambda_{\mathrm{spat-cons}}=0.05$. These hyperparameters are optimized using Bayesian optimization and kept fixed for each dataset unless otherwise specified.

Evaluation Metrics. The classification performance is evaluated using Overall Accuracy (OA), Average Accuracy (AA), and Cohen’s Kappa coefficient ($\kappa$), which are the standard metrics for HSIC. To provide a more comprehensive assessment, we additionally report the Macro-F1 score and Matthews Correlation Coefficient (MCC). Macro-F1 evaluates the harmonic mean of precision and recall averaged across all classes, whereas MCC provides a balanced measure of classification quality even in highly imbalanced scenarios.

Model calibration is assessed using the Expected Calibration Error (ECE) and multiclass Brier Score [56, 57]. ECE measures the discrepancy between prediction confidence and empirical accuracy, while the Brier Score quantifies the mean squared error between predicted probabilities and GT labels.

Lower values indicate better-calibrated probabilistic predictions. To evaluate the contribution of latent rollout reasoning, branch-wise analysis is performed by separately measuring the performance of the base classifier, trajectory classifier, and their fused prediction. In addition to OA, AA, and Kappa, we report the correction rate (CR), harm rate (HR), and net gain (NG). The CR measures the proportion of samples incorrectly classified by the base branch but correctly classified by the trajectory branch, whereas the HR measures the opposite case. Their difference, referred to as NG, quantifies the overall benefit introduced by latent rollout reasoning. Log Loss measures the quality of the predicted class probabilities by computing the average negative log-probability assigned to the GT class. Top3 accuracy measures the proportion of test samples for which the GT class appears among the three classes with the highest predicted probabilities.

Latent Rollout Contribution

To investigate the contribution of latent rollout reasoning, we analyze the performance of the individual prediction branches and their fusion. Specifically, the proposed framework produces three predictions: (i) the Base branch, which relies solely on the fused transformer representation; (ii) the Trajectory branch, which performs classification using the planned latent state obtained from stochastic rollout reasoning; and (iii) the Final prediction, which combines both branches through logit fusion. This analysis provides insight into whether latent rollout modeling contributes complementary information beyond static spatial-spectral encoding. Table I reports the branch-wise classification results. The performance difference between an individual branch and the final prediction is quantified as $\Delta= \mathrm{OA}_{\mathrm{Final}} - \mathrm{OA}_{\mathrm{Branch}}$, where $\mathrm{OA}_{\mathrm{Branch}}$ denotes the OA of either the Base or Trajectory branch. A positive $\Delta$ indicates that fusion improves upon the corresponding individual predictor.

Branch	OA	AA	$\kappa$	$\Delta$
DFC18
Base	0.9100	0.9198	0.8858	+0.0149
Trajectory	0.9113	0.9273	0.8976	+0.0136
Final	0.9249	0.9352	0.9021	–
Dioni
Base	0.8648	0.8835	0.8343	+0.0523
Trajectory	0.9037	0.9208	0.8837	+0.0134
Final	0.9171	0.9346	0.8979	–
Loukia
Base	0.8470	0.8596	0.8207	+0.0437
Trajectory	0.8811	0.8678	0.8613	+0.0096
Final	0.8907	0.8775	0.8709	–
HH
Base	0.9703	0.9643	0.9677	+0.0211
Trajectory	0.9810	0.9763	0.9787	+0.0104
Final	0.9914	0.9864	0.9891	–
TD
Base	0.9852	0.9859	0.9846	+0.0100
Trajectory	0.9946	0.9953	0.9939	+0.0006
Final	0.9952	0.9958	0.9945	–

Table I shows that the Trajectory branch consistently outperforms the Base branch across all evaluated datasets. This observation indicates that latent rollout reasoning captures predictive structure that is not fully represented by static transformer features alone. The improvement is particularly pronounced on Dioni and Loukia, where the trajectory branch increases OA from 86.48% to 90.37% and from 84.70% to 88.11%, respectively. These gains suggest that stochastic latent evolution is especially beneficial for scenes exhibiting large intra-class variability and complex spatial organization. The fused prediction consistently achieves the highest performance. The largest fusion gain is observed on Dioni, where the final model improves OA by 5.23 percentage points relative to the Base branch, highlighting the importance of latent rollout reasoning in challenging remote sensing scenarios. Although the gains are smaller on HH due to its already saturated performance, the fused prediction still improves upon both individual branches.

Beyond branch-level accuracy, we further evaluate the complementarity between the Base and Trajectory branches. Let $y_i$ denote the GT label, $\hat{y}_i^{(b)}$ the prediction of the Base branch, and $\hat{y}_i^{(t)}$ the prediction of the Trajectory branch. The CR measures the proportion of samples that are misclassified by the Base branch but correctly classified by the Trajectory branch

where $\mathbb{I}(\cdot)$ is the indicator function and $N_{\mathrm{test}}$ is the number of test samples. Similarly, the HR measures the proportion of samples correctly classified by the Base branch but misclassified by the Trajectory branch

The overall benefit of the trajectory branch is quantified through the NG as $\mathrm{NG} = \mathrm{CR} - \mathrm{HR}$. To estimate the maximum achievable benefit from combining both branches, we define the Oracle Accuracy $\mathrm{OA}_{\mathrm{oracle}}$

The Agreement Rate (AGR) measures the fraction of samples on which both branches produce identical predictions

Finally, the diversity between the predictive distributions of the two branches is measured using the Jensen-Shannon divergence (JSD) as $\mathrm{JSD} = \frac{1}{2}D_{\mathrm{KL}}(\mathbf{p}^{(b)} \mathbin\Vert\mathbf{m}) + \frac{1}{2}D_{\mathrm{KL}}(\mathbf{p}^{(t)} \mathbin\Vert\mathbf{m})$, $\mathbf{m} = \frac{1}{2}(\mathbf{p}^{(b)} + \mathbf{p}^{(t)})$ and the Kullback-Leibler divergence is defined as $D_{\mathrm{KL}}(\mathbf{p} \mathbin\Vert\mathbf{q}) = \sum_{k=0}^{K-1} p_k \log\left(\frac{p_k}{q_k}\right)$, where $\mathbf{p}^{(b)}$ and $\mathbf{p}^{(t)}$ denote the softmax probability vectors from the Base and Trajectory branches, respectively. A larger JSD indicates greater diversity between the Base and Trajectory predictions.

The correction analysis in Table II provides a direct assessment of the complementary behavior between the Base and Trajectory branches. In most datasets, the CR exceeds the HR, resulting in positive NG values. The largest gains are observed on Dioni and Loukia, where the trajectory branch corrects 6.41% and 6.37% of the errors made by the Base branch while introducing only 1.51% and 1.96% additional errors, respectively. This yields NGs of 4.89% and 4.41%, demonstrating substantial complementary information from the rollout branch on these datasets. In contrast, for TD, the CR is 0.11% while the HR is 0.18%, resulting in a negative NG of −0.06%. This indicates that the trajectory branch introduces slightly more additional errors than the number of Base-branch errors it corrects on this dataset.

Dataset	CR (%)	HR (%)	NG (%)	OAoracle (%)	AGR (%)	JSD
DFC18	1.76	1.25	0.51	93.67	96.31	0.0196
Dioni	6.41	1.51	4.89	92.89	90.86	0.0817
Loukia	6.37	1.96	4.41	91.08	90.74	0.0433
HH	0.25	0.17	0.08	99.28	99.54	0.0033
TD	0.11	0.18	-0.06	99.64	99.70	0.0017

The $OA_{\mathrm{oracle}}$ further reveals the potential upper bound obtainable from combining both branches, reaching 92.89% on Dioni and 93.67% on DFC18. Moreover, the JSD correlates strongly with the observed gains. Dioni exhibits the largest JSD (0.0817) and the largest NG, indicating that greater predictive diversity between branches translates into larger performance improvements after fusion. In contrast, HH exhibits very high AGR (99.54%) and very small JSD (0.0033), indicating that both branches converge to nearly identical predictions on this near-saturated dataset.

Robustness Analysis

A practical HSIC system should remain reliable under incomplete spectral observations and sensor noise. To evaluate the robustness of the proposed LWM framework, we consider two challenging perturbation scenarios: (i) random spectral band removal and (ii) additive Gaussian noise. For spectral corruption, a fraction $\rho$ of spectral bands is randomly removed and set to zero. Let $\mathbf{m} \in\{0,1\}^{C}$ be a random binary spectral mask satisfying $\Pr(m_c = 0) = \rho$. The mask is expanded across the batch and spatial dimensions as $\widetilde{\mathbf{m}} = \mathbf{1}_B \otimes\mathbf{m} \otimes\mathbf{1}_P \otimes\mathbf{1}_P \in\{0,1\}^{B \times C \times P \times P}$, where $\mathbf{1}_B$ and $\mathbf{1}_P$ denote vectors of ones, and $\otimes$ denotes replication along the corresponding dimensions. The corrupted input is then defined as $\widetilde{\mathbf{x}} = \mathbf{x} \odot\widetilde{\mathbf{m}}$.

For noise robustness, zero-mean Gaussian noise is added to the input $\widetilde{\mathbf{x}} = \mathbf{x} + \boldsymbol{\eta}$, $\boldsymbol{\eta} \sim\mathcal{N}(0,\sigma^2\mathbf{I})$, where $\boldsymbol{\eta} \in\mathbb{R}^{B \times C \times P \times P}$ has the same dimensionality as $\mathbf{x}$. The experiments are performed without retraining the model. Performance is evaluated using OA, AA, and $\kappa$ under increasing corruption levels. The results are reported in Table III.

The robustness results reveal several important characteristics of the proposed LWM framework. First, LWM exhibits strong resilience to additive Gaussian noise across all datasets. Even at $\sigma= 0.05$, the degradation in OA remains negligible. For example, on DFC18, the OA decreases from 92.49% to 92.42%, while on HH, the performance remains virtually unchanged. This behavior indicates that the learned latent representations are stable under moderate perturbations and are not overly dependent on individual spectral measurements.

Second, performance degrades more gradually under spectral band removal than would typically be expected for models relying heavily on individual wavelengths. Although accuracy decreases as the fraction of removed bands increases, the model retains strong predictive capability even when 30% of the spectral information is discarded. On Dioni and Loukia, the

Dataset	Perturbation	Level	OA	AA	$\kappa$
DFC18	Missing Bands	10%	0.9132	0.9341	0.8872
DFC18	Missing Bands	20%	0.8752	0.9130	0.8389
DFC18	Missing Bands	30%	0.7925	0.8514	0.7355
DFC18	Gaussian Noise	$\sigma= 0.01$	0.9249	0.9352	0.9020
DFC18	Gaussian Noise	$\sigma= 0.03$	0.9247	0.9351	0.9018
DFC18	Gaussian Noise	$\sigma= 0.05$	0.9242	0.9344	0.9012
Dioni	Missing Bands	10%	0.9167	0.9092	0.8967
Dioni	Missing Bands	20%	0.8928	0.8662	0.8661
Dioni	Missing Bands	30%	0.8372	0.7869	0.7937
Dioni	Gaussian Noise	$\sigma= 0.01$	0.9170	0.9346	0.8978
Dioni	Gaussian Noise	$\sigma= 0.03$	0.9168	0.9345	0.8975
Dioni	Gaussian Noise	$\sigma= 0.05$	0.9165	0.9343	0.8972
Loukia	Missing Bands	10%	0.8876	0.8666	0.8669
Loukia	Missing Bands	20%	0.8717	0.8435	0.8476
Loukia	Missing Bands	30%	0.8442	0.8012	0.8143
Loukia	Gaussian Noise	$\sigma= 0.01$	0.8903	0.8770	0.8704
Loukia	Gaussian Noise	$\sigma= 0.03$	0.8911	0.8782	0.8714
Loukia	Gaussian Noise	$\sigma= 0.05$	0.8905	0.8779	0.8707
HH	Missing Bands	10%	0.9917	0.9860	0.9895
HH	Missing Bands	20%	0.9887	0.9792	0.9857
HH	Missing Bands	30%	0.9686	0.9476	0.9603
HH	Gaussian Noise	$\sigma= 0.01$	0.9914	0.9865	0.9892
HH	Gaussian Noise	$\sigma= 0.03$	0.9914	0.9865	0.9892
HH	Gaussian Noise	$\sigma= 0.05$	0.9914	0.9864	0.9892
TD	Missing Bands	10%	0.9925	0.9940	0.9915
TD	Missing Bands	20%	0.9802	0.9838	0.9774
TD	Missing Bands	30%	0.9441	0.9499	0.9364
TD	Gaussian Noise	$\sigma= 0.01$	0.9952	0.9958	0.9946
TD	Gaussian Noise	$\sigma= 0.03$	0.9952	0.9958	0.9945
TD	Gaussian Noise	$\sigma= 0.05$	0.9952	0.9958	0.9945

OA remains above 83%, while HH retains an OA of 96.86%. These results suggest that the latent rollout mechanism learns distributed spatial-spectral representations that remain informative even when substantial portions of the spectrum are unavailable.

Third, the robustness trends are consistent with the design of LWM. The masked reconstruction objectives encourage recovery of missing information from surrounding spectral and spatial context, while the spectral and spatial dynamics constraints enforce structured relationships between latent representations. In addition, stochastic latent rollout reasoning promotes representation consistency across multiple plausible latent states rather than relying on a single deterministic embedding. Consequently, the model becomes less sensitive to localized corruption and measurement noise.

Attention-based Interpretability Analysis

To better understand the decision-making behavior, we analyze the spatial distribution of the learned attention responses. While quantitative metrics evaluate classification performance, they provide limited insight into how the model utilizes spatial-spectral information during inference. To address this limitation, the learned attention responses are projected back to the image plane and visualized as attention heatmaps and attention overlays on the corresponding GT scenes. These visualizations (Figs. 2–6) are further compared with classification error maps to examine the relationship between attention localization and prediction quality. Such analysis provides qualitative evidence regarding whether the latent rollout mechanism guides the model toward discriminative and spatially coherent scene representations.

The attention visualizations reveal several consistent patterns. First, high-attention responses are concentrated within semantically homogeneous regions rather than being uniformly distributed across the scene. In DFC18 (Figure 2), the learned attention focuses on dominant urban structures, including transportation corridors, vegetation regions, and built-up areas. For the Dioni and Loukia datasets (Figure 3 and Figure 4), attention responses remain spatially coherent over large land-cover regions, indicating that the model successfully captures long-range contextual relationships despite substantial intra-class variability. Similarly, for the agricultural HH dataset (Figure 5), the attention maps exhibit highly compact and class-consistent patterns that align closely with crop distributions.

Second, the attention overlays demonstrate strong correspondence between the learned responses and semantically meaningful structures in the GT scenes. Regions receiving close attention generally coincide with dominant land-cover patterns and object boundaries that are critical for discrimination. This behavior suggests that the latent rollout mechanism encourages the model to focus on informative spatial-spectral structures while suppressing irrelevant background responses.

Third, the error maps show that the majority of misclassifications occur near class boundaries and transition regions where neighboring categories exhibit similar spectral characteristics. In contrast, regions associated with strong attention responses are predominantly classified correctly.

Ablation Study

To quantify the contribution of each component, we progressively construct the proposed LWM framework from a baseline Spatial-Spectral Transformer (SST). Specifically, five configurations are evaluated: (i) SST only, which contains only the spatial-spectral transformer backbone; (ii) SST+Mask Recon, which additionally introduces masked spectral and spatial reconstruction objectives; (iii) SST+World, which incorporates the stochastic LWM module and rollout reasoning; (iv) LWM w/o Physical Priors, which further includes latent rollout planning and trajectory supervision but excludes SAM, spectral smoothness, and spatial consistency constraints; and (v) the full LWM model. The relative improvement of the full LWM model over the SST baseline is computed as Relative Gain(%) $= \frac{\mathcal{M}_{\mathrm{LWM}}-\mathcal{M}_{\mathrm{SST}}}{\mathcal{M}_{\mathrm{SST}}} \times100$, where $\mathcal{M}$ denotes the evaluation metric of interest. For metrics where lower values are better, such as Log Loss, negative values indicate performance improvement.

From Table IV, several important observations can be made. First, masked spatial-spectral reconstruction consistently improves performance over the SST baseline. Although the gains in OA are relatively modest, reconstruction supervision produces noticeable improvements in AA, $\kappa$, MCC, and calibration metrics. This indicates that recovering masked spectral bands and spatial pixels encourages the encoder to preserve richer spatial-spectral information, resulting in more robust latent representations under limited supervision.

Dataset	Setting	Loss↓	OA↑	AA↑	$\kappa$↑	MCC↑	Log Loss↓	Top3↑	Macro-F1↑	Brier↓	ECE↓
DFC18	SST only	0.5794	0.8992	0.9083	0.8683	0.8689	0.3925	0.9830	0.8919	0.1447	0.0385
DFC18	SST+Mask Recon	0.6267	0.9076	0.9190	0.8796	0.8802	0.3390	0.9841	0.8994	0.1630	0.0493
DFC18	SST+World	0.5196	0.9175	0.9324	0.8925	0.8931	0.2832	0.9878	0.9030	0.1260	0.0280
DFC18	LWM w/o Physical	0.4601	0.9225	0.9325	0.8991	0.8996	0.2942	0.9877	0.9082	0.1234	0.0296
DFC18	Full LWM	0.4100	0.9249	0.9352	0.9021	0.9024	0.2783	0.9899	0.9226	0.1192	0.0283
DFC18	Relative Gain (%)	-29.23	+2.86	+2.96	+3.89	+3.86	-29.10	+0.70	+3.44	-17.62	-26.49
Dioni	SST only	0.7405	0.8388	0.8696	0.8035	0.8069	0.4981	0.9829	0.8061	0.2341	0.0436
Dioni	SST+Mask Recon	0.7712	0.8688	0.8888	0.8395	0.8408	0.4318	0.9847	0.8436	0.2017	0.0577
Dioni	SST+World	0.6116	0.8967	0.9166	0.8731	0.8748	0.3362	0.9922	0.8625	0.1691	0.0356
Dioni	LWM w/o Physical	0.5774	0.8920	0.8915	0.8677	0.8689	0.3702	0.9850	0.9014	0.1626	0.0183
Dioni	Full LWM	0.5519	0.9171	0.9346	0.8979	0.8986	0.2894	0.9943	0.9146	0.1332	0.0080
Dioni	Relative Gain (%)	-25.47	+9.34	+7.48	+11.75	+11.36	-41.90	+1.16	+13.46	-43.10	-81.65
Loukia	SST only	0.7106	0.8541	0.8858	0.8293	0.8308	0.4360	0.9899	0.8426	0.2229	0.0415
Loukia	SST+Mask Recon	0.7201	0.8557	0.8725	0.8307	0.8315	0.4579	0.9856	0.8400	0.2038	0.0411
Loukia	SST+World	0.7296	0.8695	0.8626	0.8465	0.8470	0.4438	0.9820	0.8430	0.2021	0.0516
Loukia	LWM w/o Physical	0.6806	0.8697	0.8718	0.8471	0.8479	0.4201	0.9891	0.8446	0.2208	0.0374
Loukia	Full LWM	0.6112	0.8907	0.8775	0.8709	0.8713	0.3699	0.9902	0.8737	0.1723	0.0492
Loukia	Relative Gain (%)	-14.00	+4.29	-0.94	+5.02	+4.88	-15.16	+0.03	+3.69	-22.70	+18.55
HH	SST only	0.2452	0.9874	0.9805	0.9841	0.9841	0.0553	0.9992	0.9649	0.0252	0.0114
HH	SST+Mask Recon	0.2694	0.9906	0.9841	0.9882	0.9882	0.0484	0.9992	0.9725	0.0197	0.0119
HH	SST+World	0.1608	0.9840	0.9836	0.9799	0.9799	0.0650	0.9989	0.9739	0.0198	0.0118
HH	LWM w/o Physical	0.1759	0.9878	0.9812	0.9846	0.9846	0.0556	0.9992	0.9774	0.0152	0.0130
HH	Full LWM	0.1496	0.9914	0.9864	0.9891	0.9891	0.0486	0.9991	0.9826	0.0144	0.0131
HH	Relative Gain (%)	-38.99	+0.41	+0.60	+0.51	+0.51	-12.12	-0.01	+1.83	-42.86	+14.71
TD	SST only	0.1263	0.9893	0.9930	0.9878	0.9879	0.0652	0.9995	0.9848	0.0096	0.0148
TD	SST+Mask Recon	0.1564	0.9904	0.9926	0.9891	0.9891	0.0533	0.9996	0.9927	0.0128	0.0147
TD	SST+World	0.2661	0.9918	0.9946	0.9907	0.9907	0.0410	0.9999	0.9911	0.0191	0.0112
TD	LWM w/o Physical	0.2091	0.9942	0.9959	0.9934	0.9934	0.0364	0.9995	0.9913	0.0162	0.0141
TD	Full LWM	0.1166	0.9952	0.9958	0.9945	0.9945	0.0324	0.9999	0.9947	0.0079	0.0153
TD	Relative Gain (%)	-7.68	+0.60	+0.28	+0.68	+0.67	-50.31	+0.04	+1.01	-17.71	+3.38

Second, introducing the stochastic LWM yields substantially larger improvements than reconstruction alone. For example, on Dioni, OA increases from 83.88% to 89.67%, while $\kappa$ improves from 0.8035 to 0.8731. Similar trends are observed on DFC18 and Loukia. These results demonstrate that latent rollout reasoning provides information beyond conventional transformer-based feature aggregation.

Third, trajectory planning and trajectory-aware supervision further improve performance. Comparing SST+World with LWM w/o physical priors shows consistent gains in OA, $\kappa$, MCC, and Macro-F1 on most datasets. This observation

validates the effectiveness of aggregating multiple rollout states through learned importance weighting rather than relying solely on the initial latent representation. The results suggest that trajectory-aware aggregation enables the classifier to exploit complementary information distributed across the latent rollout trajectory.

Finally, the full LWM model achieves the best overall performance across almost all datasets. The inclusion of physically motivated priors, including spectral angle alignment, spectral smoothness, and spatial consistency regularization, consistently improves classification accuracy while simultaneously reducing probabilistic uncertainty. In particular, Dioni exhibits the largest gains, where the full model achieves relative improvements of 9.34%, 7.48%, 11.75%, and 11.36% in OA, AA, $\kappa$, and MCC, respectively, compared with the SST baseline. Furthermore, the Brier score is reduced by 43.10% and ECE by 81.65%, indicating substantially improved prediction calibration. Similar trends are observed on DFC18, where calibration improves by 17.62% and 26.49% in terms of Brier score and ECE, respectively. The HH dataset exhibits comparatively smaller improvements because the baseline transformer already achieves near-saturated performance.

Comparison with State-of-the-Art

To evaluate the effectiveness of the proposed LWM framework, we compare it against several recent state-of-the-art transformer-based HSIC methods, including Spectral-Guided Multiscale Feature-Aware Transformer (SMFAT) [58], WaveFormer (WF) [59], Differential Spatial-Spectral Transformer (DF) [60], Mid-Range Convolutional Modulated Transformer Network (MCMTN) [61], Multiscale spatial-spectral Supertoken Aggregation Transformer (MSST) [62], Spatial-Spectral Coupling and Interactive Fusion Transformer Network (S2CIFT) [63], and Spatial-Spectral CNN-Transformer Network (S2CAT) [64]. All methods are evaluated under identical experimental settings. Quantitative comparisons are reported in Table V, while qualitative error-map comparisons are presented in Figure 7-Figure 11.

Dataset	Metric	LWM	WF	DF	MCMTN	SMFAT	MSST	S2CIFT	S2CAT
Training Samples (%)		10	50	50	50	50	50	50	50
DFC18	OA	92.4997±0.1096	91.4101±0.5559	89.6109±0.0707	90.4070±0.1673	91.9506±0.1864	87.0541±0.3509	91.6866±0.0159	91.5293±0.2206
	AA	93.5275±0.0126	88.9789±0.1071	81.3057±0.3993	88.0639±0.1983	85.6948±0.1386	78.4143±0.7389	90.9763±0.3494	88.9707±0.9180
	$\kappa$	90.2191±0.0231	89.9807±0.7378	86.2516±0.1730	87.2409±0.1456	89.0046±0.1064	82.7725±0.4699	89.6388±0.0193	89.1599±0.3376
Loukia	OA	89.2107±0.1094	88.3205±0.2095	88.0258±0.9530	80.4684±0.2983	82.7826±0.1956	76.2885±1.2358	87.3655±0.1885	87.2983±0.3037
	AA	87.9875±0.1453	87.2253±0.3477	83.2571±2.2707	80.0375±0.2091	82.5046±0.2193	67.3756±0.2072	86.2443±0.1001	86.1633±0.5295
	$\kappa$	87.9916±0.1312	87.0363±0.2558	85.7381±1.1069	76.5459±0.1659	79.4113±0.2510	71.8070±1.4110	86.7121±0.2254	87.0125±0.3645
Training Samples (%)		2	2	2	2	2	2	2	2
Dioni	OA	91.7132±0.1453	90.8979±0.4669	90.7228±0.0196	88.9929±0.2319	83.2195±0.1639	89.0911±1.4596	90.9907±0.8240	90.9895±0.1569
	AA	93.4698±0.1652	89.7330±2.0580	85.1594±0.7765	76.1774±0.1982	70.9206±0.2197	75.2173±6.3917	86.8359±1.9078	89.5406±3.6479
	$\kappa$	89.7989±0.1452	88.6610±0.5898	88.9784±0.0236	86.2635±0.1596	78.9040±0.1932	86.3184±1.9485	89.5262±1.0279	88.7806±0.2026
HH	OA	99.1497±0.0914	96.9677±0.1979	96.8488±0.0625	98.1871±0.1934	98.8740±0.1864	97.1105±0.2391	97.8370±0.3560	97.0178±0.0101
	AA	98.6489±0.0912	92.5037±1.0111	91.8830±0.0323	95.3153±0.1684	97.5921±0.1274	90.6244±2.2489	94.3517±1.8307	93.0144±0.4234
	$\kappa$	98.9198±0.0731	96.1681±0.2495	96.0204±0.0793	97.7084±0.1337	98.5767±0.1425	96.3459±0.3048	97.2674±0.4510	96.2311±0.0129
TD	OA	99.5243±0.0176	96.8847±0.1554	96.9727±1.2515	98.1632±0.2984	98.7562±0.1973	97.3469±0.0576	97.8779±0.6498	96.7219±0.0587
	AA	99.5853±0.0912	93.8958±0.8300	93.5081±4.6320	98.0945±0.2532	98.6965±0.1673	93.8685±0.7378	95.2817±1.0258	93.5659±1.1773
	$\kappa$	99.4543±0.0154	96.4494±0.1780	96.5505±1.4312	97.9096±0.1903	98.5843±0.1437	96.9757±0.0664	97.5839±0.7441	96.2641±0.0634

Quantitative Results Discussion. The quantitative results in Table V demonstrate that the proposed LWM framework consistently achieves the strongest overall performance. For the Dioni dataset, which contains substantial intra-class variability and complex spatial structures, LWM achieves an OA of 91.71%, outperforming the strongest competing method by

It is important to note that the compared methods were originally designed and evaluated under different training protocols and dataset complexities. Preliminary experiments showed that several baseline methods experienced severe performance degradation and unstable convergence when trained using the same extremely limited number of labeled samples employed by LWM on the more challenging benchmark datasets. Despite using substantially fewer labeled samples, LWM consistently achieves competitive performance across all datasets.

approximately 0.7-0.8%. More importantly, the improvement in AA is considerably larger, demonstrating that the proposed framework improves classification performance not only for dominant categories but also for minority classes. A similar trend is observed on the Loukia dataset. Although several competing methods achieve competitive OA values, LWM consistently provides superior OA, AA, and $\kappa$ scores.

For the HH and TD datasets, which already exhibit extremely high baseline accuracies, LWM still achieves the best overall performance. On HH, the proposed framework exceeds all competing methods, demonstrating that latent trajectory reasoning remains beneficial even in near-saturated classification scenarios. On TD, LWM achieves an OA of 99.52%, indicating exceptional discrimination capability under highly structured coastal and urban environments. The relatively small standard deviations further demonstrate the stability and robustness of the proposed learning framework.

Qualitative Results Discussion. The qualitative comparisons further support the quantitative findings. Across all evaluated datasets (Figure 7-Figure 11), the error maps generated by LWM contain noticeably fewer misclassified pixels than competing methods. The improvement is particularly evident near class boundaries, heterogeneous regions, and spectrally ambiguous areas, where traditional attention-based models frequently produce fragmented predictions. In contrast, LWM generates smoother and more spatially coherent classification maps while preserving fine structural details. Unlike conventional transformer architectures that perform classification using a single static representation, LWM augments representation learning with stochastic latent rollout reasoning. The latent trajectory provides multiple plausible future representations that are aggregated through the rollout planner, enabling the classifier to exploit structured predictive information before making a decision.

Limitations and Failure Analysis

Although the proposed LWM framework demonstrates consistent improvements across multiple hyperspectral benchmarks, it is important to analyze its limitations and identify scenarios where latent rollout reasoning may provide limited benefit. Understanding such failure modes is essential for assessing the practical applicability of the framework and for guiding future research on latent predictive modeling for HSI analysis. In this section, we discuss computational considerations, modeling limitations, and representative failure cases observed during experimentation.

Modeling Limitations

The proposed framework introduces additional computational complexity compared with conventional transformer-based classifiers. In addition to spectral, spatial, and fusion transformers, LWM performs multi-step latent rollouts and trajectory aggregation. Although the rollout horizon used in this work is relatively short ($T = 3$), the computational cost grows approximately linearly with the rollout length. Consequently, very long rollout horizons may provide diminishing returns while increasing memory consumption and inference latency.

Another limitation stems from the absence of explicit physical dynamics. Unlike world models used in temporal environments, the latent transitions learned by LWM do not correspond to true temporal evolution. Instead, they represent structured transformations in representation space. While this formulation is effective for classification, the learned latent trajectories should not be interpreted as physically meaningful spectral or environmental changes. Consequently, the rollout process primarily acts as a predictive regularization mechanism rather than a physical simulation model.

Furthermore, the current implementation assumes a fixed patch size and local spatial context. Extremely large-scale structures or long-range spatial dependencies extending beyond the extracted patch may not be fully represented. Although transformer attention partially alleviates this limitation, global scene-level reasoning remains constrained by the local patch formulation adopted in most HSIC benchmarks.

Finally, the effectiveness of latent rollout reasoning appears to depend on dataset complexity. As observed in the branch-wise analysis, highly saturated datasets such as HH and TD exhibit extremely high agreement between the Base and Trajectory branches. In such cases, the predictive diversity between branches becomes very small, reducing the potential gains obtainable through latent rollout reasoning.

Failure Analysis

Detailed inspection of classification maps, attention visualizations, and error maps reveals several recurring failure patterns. First, most misclassifications occur near class boundaries where neighboring categories exhibit similar spectral signatures. Such boundary regions often contain mixed pixels whose spectral responses represent combinations of multiple materials. Although latent rollout reasoning improves contextual understanding, the ambiguity inherent in mixed pixels remains challenging.

Second, errors are frequently observed among semantically related classes that possess strong spectral overlap. In agricultural scenes, crops at similar growth stages may produce nearly indistinguishable spectral responses, while in urban environments different materials can exhibit highly correlated reflectance patterns. These cases provide limited discriminative evidence even after latent trajectory aggregation.

Third, the latent rollout mechanism occasionally propagates uncertainty from an incorrect initial representation. Because rollout generation begins from the environment state $z_0$, severe errors in the initial representation may lead to trajectories that remain consistently biased. This behavior is most likely to occur in highly corrupted regions, rare classes with limited training samples, or samples located near complex class transitions.

The branch correction analysis provides additional insight into these failure modes. Datasets such as Dioni and Loukia exhibit relatively large JSD and positive NG values, indicating that the trajectory branch successfully corrects a substantial number of Base-branch errors. In contrast, HH and TD show extremely high agreement rates and very small divergence values, suggesting that both branches converge to nearly identical predictions. In such scenarios, the latent rollout branch contributes limited corrective information because the underlying classification problem is already nearly solved by the static representation. Despite these limitations, the observed failure cases are largely concentrated in ambiguous boundary regions, spectrally overlapping categories, and highly challenging samples.

Conclusions and Future Research Directions

This paper introduced Stochastic Latent World Modeling (LWM), a trajectory-aware representation learning framework for hyperspectral image classification (HSIC). Unlike conventional approaches that rely on static spatial-spectral representations, the proposed framework models structured latent evolution through stochastic rollout reasoning in representation space. By constructing an environment state from fused spatial-spectral transformer features and propagating it through a stochastic latent transition process, LWM enables the classifier to exploit predictive latent structure beyond the information available in a single forward-pass representation.

Extensive experiments conducted on diverse hyperspectral benchmarks demonstrated the effectiveness of the proposed approach. The results showed that LWM consistently achieves superior classification performance compared with recent state-of-the-art methods while also improving calibration quality, robustness to spectral corruption and noise. Branch-wise analysis further confirmed that latent rollout reasoning contributes complementary predictive information beyond static transformer features. In addition, ablation studies verified the importance of each component of the framework.

Several research directions remain open for further investigation, as already discussed in Section IX. First, the current rollout mechanism operates at a single spatial scale; extending LWM to hierarchical and multi-scale representations may further improve the modeling of complex scene structures. Second, incorporating adaptive rollout horizons and uncertainty-aware planning strategies could allow the model to dynamically adjust its latent reasoning process according to scene complexity. Third, integrating LWM with semi-supervised, self-supervised, and foundation-model paradigms may improve scalability under extremely limited supervision. Finally, extending the proposed framework to hyperspectral change detection, multimodal remote sensing fusion, domain adaptation, and spatiotemporal Earth observation analysis represents an exciting direction for future research.

References