To improve the classification accuracies of hyperspectral images (HSIs), especially when using a limited number of training samples, a novel superpixel-based multiple statistical feature extraction (SPMSFE) method is proposed in this article. For each dimension-reduced pixel obtained by maximum noise fraction (MNF), the most similar superpixel-based neighbors of different sizes are first identified based on the spatial structures of the HSIs to exploit contextual spatial information accurately. Then, multiple statistical features, including the mean, covariance descriptor, and the Gaussian feature, are extracted for the set of superpixel-based neighbors to fully explore the spatial geometry information, tight correlations between different spectral bands, and spatial–spectral variations from different perspectives, respectively. In addition, these three statistical features of the pixels share the same size and can be utilized for uniform classification without any dimensionality obstacles even though the sizes of the superpixel-based neighbors for different pixels may be different. Next, we construct multiple kernels to map these multiple statistical features in the Euclidean and Riemannian manifold spaces to a uniform Hilbert space and embed them into a multitask kernelized sparse representation classification (MTKSRC) model. The constructed MTKSRC model provides a natural method to effectively fuse the multiple statistical features for excellent classification performance and robustness, especially when using limited numbers of training samples. The experimental results for three widely used HSI data sets demonstrate that the classification accuracy of the proposed SPMSFE method outperforms several latest and state-of-the-art classification methods by a large margin.