SRSSF: Stability-Reliability Spectral-Spatial Fusion for Hyperspectral Image Classification

1. Introduction

The classification of hyperspectral images [1] remains a challenge in remote sensing. It helps in disaster tracking and management, urban planning, mining, defence and surveillance, and more [2]. Each pixel is a vector of hundreds of values from emitted, transmitted, absorbed, or reflected electromagnetic waves off different surfaces or objects [3], which vary by material type and span different ranges: visible, near-infrared, shortwave-infrared, and thermal-infrared. It contains both spatial and spectral information that helps to classify objects and other information properly, impossible with an RGB image due to similar visual identity but different material composition. This large spatial variability of spectral signatures, paucity of data, atmospheric and sensor noise, redundant data, and high dimensionality make it difficult to extract and process meaningful information [4]. Different solutions exist that can be grouped into six different categories based on their functionality: band selection, dimensionality reduction, spectral-spatial feature extraction, pixel-wise classification, hybrid ensemble methods, and deep learning approaches. Pixel-wise classifiers like support vector machine (SVM) [5], K-nearest neighbor [6], etc., are the earliest naive classifiers used for a wide range of supervised classification tasks. However, they suffer from the curse of high dimensionality and computational cost. So, dimensionality reduction algorithms, e.g., principal component analysis (PCA) [7], linear discriminant analysis (LDA) [8], or independent component analysis (ICA) [9], and band selection algorithms like Mutual Information (MI) [10], Sequential Forward Floating Selection (SFFS) [5], L1-regularized embedded (LASSO) [11], and correlation pruning [12] are used along with pixel-wise classifiers to improve results. However, these results are unsatisfactory, as it completely ignores the spatial information - shape, structure, etc.- has high spectral variability and redundancy. Therefore, spectral-spatial classifiers such as extended morphological profile (EMP) and Support Vector Machine (SVM) [13], Gabor and SVM [14], extended attribute profile (EAP) and SVM [15], etc., have been developed to improve classification accuracy. However, these methods require large amounts of data and significant computational resources to achieve robust performance and scalability. More recently, hybrid approaches have emerged, which can be broadly divided into classical Machine learning (ML) ensembles, e.g., SVM + Random Forest (RF) [16], and machine learning + deep learning hybrids, e.g., EMP features fused with convolution neural network (CNN), EMP + CNN [17]. In parallel, purely deep learning pipelines such as 3D CNNs [18] and transformers [19] have been explored. Deep learning-based approaches are data-hungry with high computational cost; extensive GPU requirements, less interpretability, and reproducibility issues make them a poor choice. Our SRSSF framework belongs to the class of classical machine learning hybrid ensembles, distinct from ML+DL hybrids and pure deep learning pipelines. By combining stability-selected spectral experts with spatial and joint experts under reliability weighting, SRSSF achieves interpretable and reproducible gains without relying on deep neural architectures. Our main contributions:

• We propose SRSSF, a reliability-weighted mixture-of-experts combining a stability-selected spectral expert with full-spectrum spatial and joint experts, plus a compact-input ablation (SRSSF-SelectedInput) to isolate band-restriction effects.
• We design a shared validation-only band selection that fixes a single spectral budget $k^{★}$ per dataset/seed, applied consistently across models and ablations, avoiding test-set influence.
• We introduce a strict footprint-disjoint protocol for patch-based HSI classification, enforcing zero overlap between train/val/test receptive fields and using automated geometry-based seed acceptance.
• We pair benchmarks with interpretability analysis, including SHAP feature attribution, band-stability visualization, and expert-reliability reporting to reveal which bands and experts drive decisions.

2. Related Work

Pixel-wise classifier. Pixel-wise classifiers assign labels independently per pixel, ignoring spatial context. They are computationally light but often inaccurate. SVM [5] is widely used due to low data requirements and resistance to overfitting. It maps overlapping data into higher dimensions for separability and scales from binary to multi-class. However, in high-dimensional spectral spaces, sparsity and redundancy cause overfitting, poor weighting of bands, and heavy parameter search costs. KNN [6] classifies by local neighborhoods without distributional assumptions, performing well with small, irregular datasets, but suffers in high dimensions where distances lose meaning and generalization is weak. RF [20] ensembles decision trees over random feature subsets, reducing noise sensitivity, yet struggles in heterogeneous regions with overlapping spectra.

Dimensionality reduction and band selection. Pixel-wise classifiers suffer from high-dimensionality and redundant bands, so dimensionality reduction (PCA, LDA, ICA, KPCA [7,8,9,21]) and band selection (MI, correlation pruning [12], SFFS, LASSO) are widely used. PCA is unsupervised, reduces redundancy/noise, but is weak in discrimination. LDA is supervised, maximizes class separability, yet struggles with non-linear boundaries. ICA extracts independent components for spectral unmixing but is noise-sensitive and computationally heavy. KPCA extends PCA to non-linear manifolds but depends critically on kernel choice and is expensive. MI [10] captures non-linear dependencies for feature selection but requires reliable probability estimates. SFFS [5] dynamically adds/removes features, flexible but prone to overfitting and costly. LASSO [11] embeds selection during training, yielding sparse solutions well-suited to high-dimensional HSI, though sensitive to regularization and may discard correlated useful bands.

Spectral-spatial classifier. Spectral classifiers (e.g., SVM with PCA [7]) reduce dimensionality and denoise but ignore spatial context, limiting accuracy. EMP+SVM [13] adds morphological filters for spatial coherence, simple and interpretable but sensitive to structuring elements, heavy for large images, and prone to overfitting. EAP+SVM [15] generalizes EMP with attribute filters (area, deviation, moment), offering more flexibility and typically higher accuracy. Gabor+SVM [14] incorporates texture and orientation cues, effective for urban/man-made structures, but requires careful parameter tuning and is less suited to natural scenes with irregular textures.

Hybrid approaches. Spectral–spatial methods improved accuracy but lacked robustness and scalability. Hybrid approaches emerged to stabilize performance (ensembles), reduce redundancy (band selection), balance dimensionality (PCA/LDA), and exploit new paradigms (ML, DL, self-supervised, RL). Spectral–spatial gave context; hybrids added integration and resilience. Traditional ML hybrids (e.g., SVM+RF [16]) are more resilient to overfitting, while ML+DL hybrids (e.g., EMP+CNN [17]) learn local features and their interactions better. These improve accuracy under limited labels, scale well, but risk redundancy, complex pipelines, and tuning challenges.

Deep Learning. Deep learning supersedes hybrids in accuracy under abundance of labelled data, but hybrids remain relevant for interpretability and robustness under limited data. It is end-to-end learning; it requires no handcrafted features and achieves state-of-the-art accuracy on benchmark datasets. CNN [18] and transformer-based [19] architectures are typically used.

We acknowledge both hybrid categories and higher variants of deep learning pipelines, but our work focuses on classical ML ensembles. We do not compare against ML+DL hybrids or higher variants of deep learning baselines (2D/3D CNNs), as our emphasis is on reproducible, interpretable traditional methods. Higher variants of DL (2D/3D CNNs) depend on spatial overlap, making them incompatible with our footprint-disjoint evaluation.

3. Methodology

This section describes the spatially disjoint evaluation protocol used to prevent patch-level leakage, the classical and deep-learning baselines, the proposed Stability-Reliability Spectral-Spatial Fusion (SRSSF) architecture (Figure 1), the shared band-count selection procedure, the full algorithmic pipeline, and the interpretability methodology used throughout this study.

3.1. Spatially Disjoint, Leakage-Free Split Protocol

A major concern in HSI is that neighboring pixels are highly correlated and share similar spectra. When patch-based descriptors are used, these pixels overlap, leading to biased accuracy since the model has effectively seen similar data. To prevent such leakage, we adopt a strict footprint-disjoint sampling strategy at the level of patch footprints, rather than relying only on center pixels.

3.1.1. Footprint Locking

For a candidate center pixel $(r,c)$ and patch size p, definition of its footprint as the square neighborhood: Eq.(1)

$$\mathcal{F}(r,c)=\left\{(i,j):r-⌊\frac{p}{2}⌋\le i\le r+⌊\frac{p}{2}⌋,c-⌊\frac{p}{2}⌋\le j\le c+⌊\frac{p}{2}⌋\right\}\cap [0,H)\times [0,W)\}.$$

(1)

Boolean masks, $M_{train}$ and $M_{val}$, of size $H\times W$ track pixels already locked by train or validation footprint, respectively. A center $(r,c)$ is footprint-free with respect to a mask M if $M\left[\mathcal{F}\right(r,c\left)\right]$ contains no True entries.

Training and validation centers must be footprint-free relative to $M_{train}$ and $M_{val}$, ensuring no overlap between patches within or across these sets. This stricter no-overlap rule minimizes spatial redundancy in development data used for stability scoring, expert models, and fusion weights. Only test patches may overlap (Section 3.1.3), since they are excluded from fitting, band selection, and fusion-weight estimation; overlap confined to the test set does not constitute leakage.

3.1.2. Per-Class Training and Validation Selection

For each class $k\in \{1,\dots ,C\}$ and random seed S, class- and seed-specific random states are used to shuffle the coordinates of all labelled pixels of class k. Fixed fractions of $f_{tr}$ and $f_{val}$, per dataset, are used to compute a target number of training centers $n_{tr}^k$ and validation centers $n_{val}^k$ (Table 1):

$$n_{tr}^k=max\left(1,⌈|C_k|·f_{tr}⌋\right),n_{val}^k=max\left(1,⌈|C_k|·f_{val}⌋\right)$$

(2)

where $C_k$ is the set of labelled pixels of class k; validation centers are reserved first. Training centers are fetched from the remaining candidates, under the strict constraint that each of their footprints is simultaneously free of $M_{val}$ and $M_{train}$. Once accepted, a center’s footprint is locked into its corresponding mask. This guarantees that no training footprint intersects any validation footprint.

3.1.3. Strict Test Selection

A labelled pixel $(r,c)$ is assigned to the test set only if its entire footprint $\mathcal{F}(r,c)$ is free of both $M_{\mathrm{train}}$ and $M_{\mathrm{val}}$:

$$(r,c)\in \mathcal{T}\iff M_{train}\left[\mathcal{F}(r,c)\right]=⌀\wedge M_{val}\left[\mathcal{F}(r,c)\right]=⌀.$$

(3)

It is strictly stronger than selecting center pixels after excluding train and val center pixels, which would have allowed a test patch to share up to $p^2-1$ pixels with a training or validation patch. Test patches are allowed to overlap one another, since test samples never participate in parameter fitting, band selection, or fusion-weight estimation; overlap confined to the test partition does not constitute leakage.

3.1.4. Adaptive Density Relaxation

In datasets with small or spatially fragmented classes (notably KSC [22]), strict exclusion of training and validation data may leave little or no test data. To address this, target counts for training and validation are scaled by a factor $\gamma \in \{1,0.5,0.25,\dots \}$, halved iteratively (up to six attempts) until a non-empty test region is obtained. If no valid outcome is found, a relaxed fallback (center-pixel exclusion only) is applied, and the split is explicitly flagged as non-strict. Such relaxed splits are excluded from the accepted-seed pool described in Section 3.2.

3.1.5. Leakage Verification

For every candidate split, three binary footprint masks are rasterized from the accepted train, validation, and test centers, and pairwise pixel-wise intersections are computed: Eq.(4).

$${\Delta }_{tr-val}=\left|M_{train}\wedge M_{val}\right|,{\Delta }_{tr-te}=\left|M_{train}\wedge M_{test}\right|,{\Delta }_{val-te}=\left|M_{val}\wedge M_{test}\right|.$$

(4)

A split is considered leakage-free if and only if ${\Delta }_{tr-val}={\Delta }_{tr-te}={\Delta }_{val-te}=0$. This check is performed independently of, and after, the split construction, so that any residual overlap is measured directly on the realized footprints rather than assumed from the sampling procedure.

3.2. Seed Filtering for Statistically Reliable Evaluation

3.2.1. Seed Selection

Because footprint locking is greedy and order-dependent, the number of centers per partition varies across random seeds. To ensure leakage-free and statistically meaningful results, we select exactly $R=10$ acceptable seeds per dataset from candidates $s=0,1,2,\dots$. A seed is accepted only if:

1.

the split is strict (non-empty test set without density relaxation);

2.

all three pairwise footprint overlaps are zero;

3.

the number of test centers meets the dataset-specific minimum ${\tau }_d$ (Table 2).

Acceptance depends solely on split geometry and ground-truth labels; no classifiers are fit, and no test predictions are used, preventing bias. The R accepted seeds are reused consistently across all experiments (main comparison, per-class analysis, few-shot, ablation, and efficiency benchmarking).

3.3. Data, Feature, and Band-Count Sharing/Caching Protocol

To ensure fair comparisons, we separate cached data (for reuse only) from shared splits/features (for methodological fairness). Normalized cubes, labels, and global EMP/EAP/Gabor [13,14,15] maps are cached once per dataset. Within each evaluation condition (dataset, seed, experiment type), identical splits, labels, raw spectra, and shared band-count $k★$ are passed to all methods. Baseline transforms (PCA, ICA, SpaBS, MVPCA [5,7,9,23,24]) and SRSSF band indices are recomputed independently, while few-shot experiments construct distinct splits per shot count.

3.4. Spatial and Spectral Feature Extraction

Three complementary spatial feature stacks are computed for each accepted split, once per dataset from the grayscale image $\overline{\mathbf{X}}={\mathrm{mean}}_b\left(\tilde{\mathbf{X}}\right)$ and cached for reuse across all train, validation, and test centers.

Extended Morphological Profile (EMP). Disk-shaped structure elements with radii $\{1,2,3\}$ are applied with grayscale morphological opening and closing, yielding a six-channel morphological stack; mean filters with window sizes $\{3,5,7\}$ are applied as a fallback when no morphology library is available.

Extended Attribute Profile (EAP) surrogate. Grayscale images are binarized with four intensity thresholds $\{0.1,0.3,0.5,0.7\}$, and a $5\times 5$ mean filter smooths each binary mask and its complement, producing an eight-channel attribute-like stack.

Gabor filter bank. Three spatial frequencies $\{0.10,0.25,0.40\}$ and four orientations $\{0,\pi /4,\pi /2,3\pi /4\}$ are used to construct a bank of Gabor kernels, via FFT-based convolution, yielding a twelve-channel texture stack.

Spectral vectors are extracted directly as ${\mathbf{x}}_i=\tilde{\mathbf{X}}(r_i,c_i,:)\in {\mathbb{R}}^B$ for a coordinate set $\mathcal{S}$, and the corresponding rows of the cached EMP, EAP, and Gabor stacks are gathered at the same coordinates.

3.5. Classical and Deep-Learning Baselines

We evaluate 15 configurations: 13 classical baselines (7 spectral-only, 5 spatial-spectral, and 1 deep-learning) together with the two proposed SRSSF variants. All models are fit and evaluated identically on each accepted split, all fit identically on each accepted split.

Spectral-only: Raw+SVM [5] uses standardized spectra with an RBF-SVM ($C=100$, $\gamma =\mathtt{scale}$). PCA+SVM [7] and FastPCA+RF [25,26] project onto 30 principal components (standard vs. randomized SVD), followed by SVM or RF (200 trees, depth 25). ICA+SVM [9] applies FastICA with whitening, falling back to PCA if unstable. KPCA+SVM [21] uses Nyström RBF with 300 landmarks plus PCA before SVM. SpaBS+SVM [5,23] selects 30 bands by variance, while MVPCA+SVM [5,24] ranks bands by loadings of 10 FastPCA components.

Spatial-spectral: EMP+SVM [13], EAP+SVM [15], and Gabor+SVM [14] concatenate spectra with spatial features. EMP+RF [13,26] applies RF (200 trees, depth 25) to EMP-augmented inputs. SVM+RF [16] is a soft-voting ensemble of RBF-SVM and RF on raw spectra.

Deep learning: Spectral1DCNN is a lightweight 1D CNN [27] trained on standardized spectra (B bands as a length-B sequence). It has three convolutional blocks (kernels $7,5,3$; channels $32,64,128$; BN, ReLU, pooling) and a two-layer head with dropout ($p=0.3$). Training uses Adam (${10}^{-3}$ lr, ${10}^{-4}$ wd) for 60 epochs, batch size 64. If PyTorch is unavailable, Raw+SVM is substituted and logged.

3.6. Proposed SRSSF Architecture

The proposed Stability-Reliability Spectral-Spatial Fusion (SRSSF) model is a mixture-of-experts architecture that combines full-spectrum spatial and joint representations with a compact, stability-selected spectral representation, fused via reliability weights derived from evaluation on the validation set.

3.6.1. Stability-Selected Band Subset

Given training spectra ${\mathbf{X}}_{tr}\in {\mathbb{R}}^{N\times B}$ and labels ${\mathbf{y}}_{tr}$, band relevance is estimated via $N_{boot}=10$ class-stratified bootstrap resamples (sampling fraction $0.80$ per class). For each bootstrap t, a one-way ANOVA F-score $F_b^{\left(t\right)}$ is computed per band b using f_classif.

$$F_b^{\left(t\right)}={\displaystyle \frac{\mathrm{between}-\mathrm{class}\mathrm{variance}\mathrm{of}\mathrm{band}b\mathrm{on}\mathrm{bootstrap}t}{\mathrm{within}-\mathrm{class}\mathrm{variance}\mathrm{of}\mathrm{band}b\mathrm{on}\mathrm{bootstrap}t}}.$$

(5)

The mean and standard deviation of $F_b^{\left(t\right)}$ across bootstraps define a stability-adjusted score with stability penalty $\lambda =0.50$. (Eq.(6))

$${\sigma }_b={\mu }_b-\lambda {\varsigma }_b,{\mu }_b={\displaystyle \frac{1}{N_{boot}}}\sum _t{F}_b^t,{\varsigma }_b=\sqrt{{\displaystyle \frac{1}{N_{boot}}}{\sum }_t{(F_b^t-{\mu }_b)}^2}$$

(6)

with stability penalty $\lambda =0.50$. Bands are ranked in descending order of ${\sigma }_b$ and greedily accepted into a subset $\mathcal{K}$ of size k subject to a redundancy constraint: a candidate band is admitted only if

$$\underset{b^{\prime }\in \mathcal{K}}{max}\left|\rho ({\mathbf{X}}_{tr}[:,b],{\mathbf{X}}_{tr}[:,b^{\prime }])\right|<{\rho }_{max},{\rho }_{max}=0.98,$$

(7)

where $\rho (·,·)$ is the Pearson correlation coefficient; If Eq.(7) is satisfied by fewer than k bands, the redundancy filter is not used in ranking to guarantee exactly k selected bands.

3.6.2. Expert Feature Blocks

Let ${\mathbf{X}}_{\mathcal{K}}$ denote the spectrum restricted to the stability-selected band subset $\mathcal{K}$. The spectral–spatial input is defined as ${\mathbf{X}}_{\mathrm{sp}}=\mathbf{X}$ for SRSSF(Ours) and ${\mathbf{X}}_{\mathrm{sp}}={\mathbf{X}}_{\mathcal{K}}$ for SRSSF-SelectedInput. The seven expert feature blocks are

$$\mathcal{Z}=\left\{\begin{array}{cc}{\mathbf{Z}}_{\mathrm{spectral}}={\mathbf{X}}_{\mathcal{K}},\hfill & {\mathbf{Z}}_{\mathrm{emp}}=[{\mathbf{X}}_{\mathrm{sp}}\parallel \mathrm{EMP}],\hfill \\ {\mathbf{Z}}_{\mathrm{eap}}=[{\mathbf{X}}_{\mathrm{sp}}\parallel \mathrm{EAP}],\hfill & {\mathbf{Z}}_{\mathrm{gabor}}=[{\mathbf{X}}_{\mathrm{sp}}\parallel \mathrm{Gabor}],\hfill \\ {\mathbf{Z}}_{\mathrm{joint}}=[{\mathbf{X}}_{\mathrm{sp}}\parallel \mathrm{EMP}\parallel \mathrm{EAP}\parallel \mathrm{Gabor}],\hfill & {\mathbf{Z}}_{\mathrm{rf}}={\mathbf{Z}}_{\mathrm{joint}},\hfill \\ \multicolumn{2}{c}{{\mathbf{Z}}_{\mathrm{FastPCA}-\mathrm{RF}}={\mathrm{PCA}}_{30}\left({\mathbf{X}}_{\mathrm{sp}}\right).}\end{array}\right.$$

(8)

Here, ∥ denotes column-wise concatenation and ${\mathrm{PCA}}_{30}$ is randomized-SVD PCA with at most 30 components, fit using the training partition and applied unchanged at inference. Thus, SRSSF(Ours) uses selected bands only for the spectral expert while retaining the full spectrum in all spatial and joint experts. In SRSSF-SelectedInput, ${\mathbf{X}}_{\mathrm{sp}}$ is replaced by ${\mathbf{X}}_{\mathcal{K}}$, so every expert operates on the same compact selected-band representation.

3.6.3. Expert Classifiers and Reliability-Weighted Fusion

Each expert block trains an independent classifier: an RBF-SVM ($C=100$) for spectral, emp, eap, gabor, and joint; a random forest (200 trees, depth 25) for rf and fastpca_rf. Each classifier produces class-probability estimates on the held-out validation partition, and per-expert validation accuracy $a_e$ is computed. Reliability weights are obtained via a softmax over validation accuracy with temperature $T=0.05$:

$$w_e={\displaystyle \frac{exp(a_e/T)}{{\sum }_{e^{\prime }}exp(a_{e^{\prime }}/T)}}$$

(9)

At inference, per-expert class-probability vectors ${\mathbf{p}}_e\left(\mathbf{x}\right)$ are combined as

$$\mathbf{p}\left(\mathbf{x}\right)={\displaystyle \frac{{\sum }_e{w}_e{\mathbf{p}}_e\left(\mathbf{x}\right)}{{\sum }_e{w}_e}}$$

(10)

and the predicted label is $\widehat{y}\left(\mathbf{x}\right)=arg{max}_c{p}_c\left(\mathbf{x}\right)$. After the validation-based fusion weights are fixed, all expert classifiers (and, where applicable, the PCA projection) are refit on the union of validation and training data, restricted during initial fit. So the final deployed model benefits from the additional validation samples without altering the reliability weights, computed from the held-out validation predictions.

3.7. Shared Band-Count Selection

Once k stability bands are selected per seed and dataset from the train and validation partition only, we use it across all methods and ablation configurations for both variants of SRSSF. For each candidate $k\in \{15,30,45,60,75\}$, the hybrid architecture (SRSSF(Ours), select_input_mode=False) is fit on the train partition and evaluated on the validation partition, yielding validation OA$\left(k\right)$. The smallest k within a tie tolerance of $0.25$ percentage points of the best observed validation OA is selected:

$$k^{★}=min\left\{k:O{A}_{val}\left(k\right)\ge \underset{k^{\prime }}{max}O{A}_{val}\left(k^{\prime }\right)-0.25\right\}$$

(11)

Compact spectral representations are favored over larger ones if it is meaningful. This guarantees that the test partition is never touched, except at final testing. All band-count sweeps and the resulting $k^{★}$ per dataset/seed are logged for reproducibility.

3.8. Ablation Configurations

To isolate the contribution of stability-based band selection and reliability-weighted fusion, six configurations are compared under identical splits and shared $k^{★}$: (i) RawSVM, an RBF-SVM on the full raw spectrum; (ii) SRSSF_SpectralOnly_WithSelection, an RBF-SVM trained only on the stability-selected bands ${\mathbf{X}}_{\mathcal{K}}$; (iii) SRSSF_FusionNoBandSelection, the full hybrid fusion architecture with the spectral expert forced to use all B bands instead of $\mathcal{K}$; (iv) SRSSF_UniformFusion, the hybrid fusion architecture with reliability weights replaced by a uniform average $w_e=1/\left|\mathcal{E}\right|$; (v) SRSSF(Ours); and (vi) SRSSF-SelectedInput.

3.9. Interpretability and Explainability

Beyond aggregate accuracy metrics, we report three complementary interpretability analyses, each computed per dataset, seed, and experiment type, and subsequently aggregated across the R accepted seeds.

3.9.1. Shapley-Value Feature Attribution (SHAP)

For every classical baseline and for the 1D-CNN, feature-level attributions are computed using Shapley additive explanations (SHAP). For a trained probability model $f\left(\mathbf{x}\right)=Pr(y\mid \mathbf{x})$, the SHAP value ${\varphi }_j\left(\mathbf{x}\right)$ of feature j for instance $\mathbf{x}$ is defined as the Shapley value of j in the cooperative game whose value function is the model’s expected output over feature subsets:

$${\varphi }_j\left(\mathbf{x}\right)=\sum _{S\subseteq \mathcal{J}\setminus \left\{j\right\}}{\displaystyle \frac{\left|S\right|!\left(\right|\mathcal{J}|-|S|-1)!}{\left|\mathcal{J}\right|!}}\left[f_{\mathbf{x}}(S\cup \left\{j\right\})-f_{\mathbf{x}}\left(S\right)\right],$$

(12)

where $\mathcal{J}$ is the full feature index set and $f_{\mathbf{x}}\left(S\right)$ is the model output when only the features in S are set to their values from $\mathbf{x}$ and all others are marginalized using a background sample. In practice, ${\varphi }_j\left(\mathbf{x}\right)$ is estimated with a model-agnostic Explainer using a background set of 20 training instances and an explanation set of up to 50 test instances per (dataset, seed, experiment, method) combination. Per-feature importance is summarized as the mean absolute SHAP value over the explanation set,

$${\overline{\varphi }}_j={\displaystyle \frac{1}{|{\mathcal{E}}_{\mathrm{shap}}|}}\sum _{\mathbf{x}\in {\mathcal{E}}_{\mathrm{shap}}}\left|{\varphi }_j\left(\mathbf{x}\right)\right|,$$

(13)

and the top-20 features ranked by ${\overline{\varphi }}_j$ are reported, along with aggregate summaries averaged over all R accepted seeds and, separately, over all evaluated methods. SHAP attribution is not computed for the 1D-CNN’s internal convolutional activations, since perturbation-based explainers are not directly applicable to raw sequence inputs; this limitation is explicitly noted in the explainability artifacts.

3.9.2. Band-Stability Diagnostics

For each SRSSF instance, the per-band mean ${\mu }_b$, standard deviation ${\varsigma }_b$, stability score ${\sigma }_b$ (Eq. (6)), and binary selection indicator $\mathbb{1}[b\in \mathcal{K}]$ are saved and visualized as stability curves with selected bands highlighted. Aggregation across R seeds averages ${\mu }_b,{\varsigma }_b,{\sigma }_b$, and selection frequency $\frac{1}{R}{\sum }_r\mathbb{1}[b\in {\mathcal{K}}^{\left(r\right)}]$, yielding a seed-robust view of consistently discriminative spectral regions.

3.9.3. Expert Reliability Diagnostics

For each SRSSF instance, the global expert reliability weight $w_e$ (Eq. (9)) and corresponding per-class weight (identical across classes under the current scheme) are saved and visualized as bar charts over the seven experts (Eq. (8)). Aggregating $w_e$ across R seeds highlights which experts (spectral, EMP, EAP, Gabor, joint, RF, FastPCA-RF) are consistently most reliable, complementing the model-agnostic SHAP attributions.

4. Experiments

4.1. Datasets

This study employs four benchmark hyperspectral datasets: Indian Pines, Pavia University, Salinas, and Kennedy Space Center(KSC) - [22]. Indian Pines: 220 bands (200 retained), 16 crop/vegetation classes, 145×145 pixels. Pavia University: 115 bands (103 retained, 0.43–0.86µm), 9 urban classes, 610×340 pixels. Salinas: 224 bands (204 retained), 16 agricultural classes (broccoli, celery, grapes, lettuce, vineyards, fallow, stubble, soils), 512×217 pixels. KSC: 224 bands (176 retained), 13 vegetation/water classes, 512×614 pixels.

Each dataset $\mathbf{X}\in {\mathbb{R}}^{H\times W\times B}$ contains $H\times W$ spatial pixels and B spectral bands, along with ground-truth label map $\mathbf{Y}\in {\{0,1,\dots ,C\}}^{H\times W}$. C is the number of land-cover classes, where class 0 denotes unlabelled background. Each cube is min-max normalized to $[0,1]$ using the global minimum and maximum of the full cube, before any experiment:

$$\tilde{\mathbf{X}}={\displaystyle \frac{\mathbf{X}-min\left(\mathbf{X}\right)}{max\left(\mathbf{X}\right)-min\left(\mathbf{X}\right)+ϵ}},$$

(14)

$ϵ={10}^{-8}$, prevents division by zero. Details of the fixed square patch size p used for spatial feature extraction and footprint locking for each dataset d are shown in Table 3.

4.2. Evaluation Metrics

Performance is evaluated on the held-out test centers using overall accuracy (OA), average accuracy (AA), Cohen’s kappa coefficient ($\kappa$) Eq.(16), and per-class accuracy, precision, recall, and F1-score. Let $\mathbf{C}\in {\mathbb{N}}^{C\times C}$ denote the test confusion matrix, where $C_{ij}$ is the number of pixels of true class i classified as class j. OA and AA are defined as

$$\mathrm{OA}={\displaystyle \frac{{\sum }_{i=1}^C{C}_{ii}}{{\sum }_{i=1}^C{\sum }_{j=1}^C{C}_{ij}}},\mathrm{AA}={\displaystyle \frac{1}{C}}{\displaystyle \sum _{i=1}^C}{\displaystyle \frac{C_{ii}}{{\sum }_{j=1}^C{C}_{ij}}}$$

(15)

$$\kappa ={\displaystyle \frac{p_o-p_e}{1-p_e}},p_o=\mathrm{OA},p_e={\displaystyle \frac{1}{N^2}}\sum _i\left(\sum _j{C}_{ij}\right)\left(\sum _j{C}_{ji}\right)$$

(16)

Cohen’s $\kappa$ Eq.(16) is computed from the same confusion matrix to account for chance agreement. Precision, recall, and F1-score are reported separately for each class in the per-class analysis. F1 mean score is calculated as:

$$\mathrm{F}1={\displaystyle \frac{1}{C}}{\displaystyle \sum _{i=1}^C}{\displaystyle \frac{2·C_{ii}}{{\sum }_{j=1}^C{C}_{ij}+{\sum }_{j=1}^C{C}_{ji}}}$$

(17)

Structural Breakdown, $C_{ii}$: The diagonal elements (correctly classified land-cover pixels for class $i.su{m}_{j=1}^C{C}_{ij}$: The row sum for class i, representing total ground truth pixels (User’s / Producer’s base reference). $su{m}_{j=1}^C{C}_{ji}$: The column sum for class i, representing all pixels assigned to that class by the model.

All results are reported as mean ± standard deviation over the ten accepted seeds; the realized number of successful runs is also reported for each dataset–method pair.

4.3. Implementation Details

All experiments are completed on a 64-core CPU with 432 GB RAM. The global NumPy seed is fixed to 42, while each stochastic split, bootstrap resampling, classifier initialization, and CNN run is controlled by its accepted seed for exact reproducibility. Threading-sensitive libraries (OpenMP, OpenBLAS, MKL, NumExpr) are pinned to a single thread to ensure deterministic behavior. Explainability artifacts, per-band stability scores and indicators, per-expert reliability weights, and SHAP-based feature-importance summaries (Section 3.9), are saved per dataset, experiment type, and seed, then aggregated across R accepted seeds. Split-quality diagnostics, strictness flags, footprint-overlap counts (Eq. (4)), and realized train/validation/test center counts are persisted for every seed to fully audit the leakage-free protocol (Section 3.1) and caching/sharing behavior (Section 3.3).

5. Results

The proposed SRSSF framework is evaluated on four benchmark hyperspectral datasets under strict footprint-disjoint splits, with comparisons to classical ML and 1D CNN baselines. Results are reported as mean ± std over $R=10$ accepted, leakage-free seeds per dataset (Section 3.2). Section(Section 5.1) gives quantitative table outputs; Section(Section 5.2) summarizes the qualitative visual results.

5.1. Quantitative Results

5.1.1. Main Comparison

Table 4 reports OA, AA, and $\kappa$ for all fifteen methods on all four datasets. SRSSF(Ours) attains the best OA, AA and $\kappa$ on every dataset (OA: $56.43\pm 3.21$ / $91.94\pm 1.16$ / $90.50\pm 1.62$ / $75.96\pm 2.99$ on Indian Pines/Pavia University/Salinas/KSC respectively [22]), with the largest AA margin over the best baseline on KSC ($66.70$ vs. $61.91$, EMP+SVM [13]) and Indian Pines ($48.85$ vs. $43.33$, Gabor+SVM [14]), indicating gains concentrated on under-represented classes. KPCA+SVM [21] collapses on KSC ($30.15$ OA), reflecting Nyström instability on sparse label geometry. 1D-CNN [27] is weakest everywhere (e.g., $40.34$ OA on KSC), consistent with limited data under strict disjoint splits. PCA+SVM [7]/ICA+SVM [9] are near-identical (ICA falls back to PCA on convergence failure). SRSSF-SelectedInput trails SRSSF(Ours) on every dataset (e.g., $-4.20$ OA on Indian Pines) yet still beats all thirteen non-SRSSF baselines on Pavia University and KSC in OA.

5.1.2. Few-Shot Sensitivity

Table(Table 5) reports OA at 5–25 samples/class. SRSSF(Ours) leads at nearly every budget on all datasets (SRSSF-SelectedInput is not considered for comparison against SRSSF(Ours)), with the largest relative margin at the lowest budget on Pavia University ($68.2$ vs. $61.8$ for FastPCA+RF [25,26], at 5 shots/class), confirming that stability-selected spectral features and reliability fusion are most valuable when labels are scarce. Ensemble/variance baselines (SVM+RF [16], SpaBS+SVM [5,23]) close or reverse the gap only at the largest budget on Indian Pines ($71.0$/$69.4$ vs. $66.3$ at 25 shots/class) and KSC ($99.1$ vs. $97.9$), while KPCA+SVM [21] and 1D-CNN [27] are unstable at low shots on KSC (e.g., KPCA+SVM: $39.8$ OA at 5 shots/class). AA and Kappa showed consistent trends with OA across all datasets, so we report OA only in few-shot tables.

5.1.3. Computational Efficiency

Efficiency analysis in Table 6 shows that while classical methods are fastest and most memory-efficient, they sacrifice accuracy. Deep learning and the proposed SRSSF framework incur higher training cost, but inference remains sub-second. The Selected-Input variant further reduces memory and inference time, making SRSSF practical for real-world hyperspectral applications.

5.2. Qualitative Results

Split maps The split maps in Figure 2 confirm leakage-free partitions across all four datasets: green train centers, blue validation centers, and red test centers, with gray marking footprint-locked exclusion zones. Among them, KSC exhibits the sparsest and most fragmented labelled regions, reflecting the reduced training and validation counts imposed by the strict footprint-disjoint protocol. (Section(Section 3.1)).

Classification maps (Figure 3) The classification panel presents the false-color composite, grayscale image, ground truth, and per-method prediction maps for Indian Pines (seed=0). While qualitative maps under strict footprint-disjoint splitting appear fragmented, they illustrate the relative coherence of SRSSF compared to classical baselines. The figure is included for transparency, showing that visual artifacts stem from label scarcity and spatial disjointness rather than methodological flaws. It is in Table 4.

Per-class F1 heatmap (Figure 4) The heatmap consolidates per-class F1 across all fifteen methods and four datasets. SRSSF attains the highest mean per-class F1 overall ($65.5\%$, vs. $62.9\%$ for the best classical baseline, EMP+SVM) and ranks among the top three methods on $63\%$ of classes, though it is not uniformly best at the individual-class level — several spatial baselines remain competitive on specific classes. On Indian Pines, rare classes with only a handful of labelled pixels (Alfalfa, Grass-pasture-mowed, Oats) yield near-zero F1 across all methods, confirming a label-scarcity ceiling rather than a method-specific weakness. On KSC, Oak/Broadleaf and Hardwood remain the hardest classes for SRSSF (F1 $15--20\%$), pointing to spectral confusion among visually similar wetland vegetation types as the dominant remaining error source.

Band-stability curves (Figure 5) Selected bands concentrate within compact, dataset-specific spectral regions rather than spreading uniformly, with the stable sub-bands’ location and spread varying across datasets. This indicates that the stability criterion adapts to each scene’s discriminative spectral structure rather than applying a fixed selection pattern.

Expert-reliability weights (Figure 6) Across datasets, SRSSF emphasizes different experts depending on scene characteristics. In KSC, FastPCA-RF dominates, indicating that spectral dimensionality reduction captures most discriminative variance. In Indian Pines, the RF expert leads, reflecting the importance of spatial–spectral fusion in heterogeneous agricultural fields, with FastPCA-RF secondary. In Pavia University, EMP receives the highest weight, consistent with the relevance of morphological profiles for urban structures, with Joint features close behind. In Salinas, FastPCA_RF receives the largest weight, while other experts contribute diffusely, indicating that spectral compression is most reliable but spatial cues provide complementary support. These results confirm that SRSSF adaptively prioritizes the most informative experts according to dataset properties.

SHAP attributions (Figure 7) The SHAP heatmap shows that classical spectral baselines and SRSSF’s spectral expert emphasize overlapping but non-identical band regions, confirming that stability selection retains genuinely discriminative spectral cues rather than an arbitrary subset of baseline features.

Taken together, section(Section 5.2) forms a coherent narrative: SRSSF gains accuracy by enforcing leakage-free partitions, emphasizing compact stable spectral regions, and adaptively weighting spatial–spectral experts. Residual errors stem from label scarcity and spectral confusion, highlighting both the strengths and remaining challenges of the framework.

6. Ablation Study

Table(Table 7) isolates band selection and reliability fusion using six configurations, on identical splits and shared $k^{★}$.

Three findings emerge. Spectral-Only+Selection, discarding all spatial context, is the weakest SRSSF-family configuration everywhere (e.g., $-12.8$ OA vs. SRSSF(Ours) on Indian Pines), showing spatial context, not spectral compaction, drives most of SRSSF’s accuracy. Fusion-No-Band-Selection performs statistically comparably to, and marginally above, SRSSF(Ours) on all four datasets. This indicates that band selection’s main benefit is dimensionality reduction of the spectral expert at negligible accuracy cost, more consequential for the compact SRSSF-SelectedInput variant. Replacing validation-derived reliability weights with uniform averaging costs $\approx 0.24$–$0.45$ OA points on every dataset, a modest but consistent margin confirming reliability weighting reliably improves over naive averaging.

7. Discussion

Implications. Implications. SRSSF yields the best OA/AA/$\kappa$, with the largest gains under label scarcity and spectral confusion, driven primarily by reliability-weighted fusion. Band selection contributes negligible additional accuracy on its own (SRSSF_FusionNoBandSelection matches or marginally exceeds SRSSF(Ours) in the ablation). Still, it delivers the framework’s efficiency and compactness gains at essentially no accuracy cost.

Comparison. Classical spatial descriptors (EMP [13]/EAP [15]/Gabor [14]) improve over raw spectra but less than SRSSF’s learned fusion, underscoring that expert combination matters. Unlike prior work, we enforce footprint-disjoint splits and report variance across seeds for stricter evaluation.

Limitations. Gains incur up to $\approx 128\times$ fit-time and $\approx 17\times$ memory (excluding 1D-CNN). Extreme scarcity (<10 pixels) remains unsolved; global weighting may under-serve classes; cross-scene generalization untested.

Efficiency. SVM [5] baselines are fastest but saturate early; ensembles add moderate cost with modest gains. 1D-CNN [27] and SRSSF are heavier, but inference is sub-second. SelectedInput reduces memory and inference cost while retaining accuracy, quantifying the accuracy–efficiency balance.

Future work. Per-class/pixel weighting, semi-supervised augmentation, cross-scene transfer, and pruning expert subsets.

8. Conclusion

We introduced SRSSF, a reliability-weighted mixture-of-experts combining stability-selected spectral and full-spectrum spatial/joint experts, evaluated under strict footprint-disjoint splits. Across four benchmarks and 15 methods, SRSSF achieved best OA/AA/$\kappa$, with largest gains on scarce/confusable scenes at bounded training-time and memory cost. Ablations show fusion drives accuracy, while reliability weighting and band selection improve efficiency and consistency. Remaining errors stem from scarcity and spectral similarity, motivating per-class weighting, semi-supervised augmentation, and cross-scene transfer.

Practicality. Despite higher fit-time and memory, SRSSF maintains sub-second inference. The SelectedInput variant further reduces overhead while retaining accuracy, demonstrating state-of-the-art performance with acceptable efficiency trade-offs.

9. Acknowledgment

The authors acknowledge Kaggle contributors Douglas Martins and Sreevalli Manda for all datasets, and Microsoft (Azure) for cloud credits for this work.