Spatial-spectral hierarchical vision permutator for hyperspectral image classification

ABSTRACT

In recent years, the convolutional neural network (CNN) has been widely applied to hyperspectral image classification because of its powerful feature capture ability. Nevertheless, the performance of most convolutional operations is limited by the fixed shape and size of the convolutional kernel, which causes CNN cannot fully extract global features. To address this issue, this article proposes a novel classification architecture named spatial-spectral hierarchical Vision Permutator (S2HViP). It contains a spectral module and a spatial module. In the spectral module, we divide the data into groups along the spectral dimension and treat each pixel within the group as a spectral token. Spectral long-range dependencies are obtained by fusing intra- and inter-group spectral correlations captured by multi-layer perceptrons (MLPs). In the spatial module, we first model spatial information via morphological methods and divide the resulting spatial feature maps into spatial tokens of uniform size. Then, the global spatial information is extracted through MLPs. Finally, the extracted spectral and spatial features are combined for classification. Particularly, the proposed MLP structure is an improved Vision Permutator, which presents a hierarchical fusion strategy aiming at generating discriminative features. Experimental results show that S2HViP can provide competitive performance compared to several state-of-the-art methods.

KEYWORDS:

• Hyperspectral image classification
• vision permutator
• feature fusion
• long-range dependencies

Introduction

Hyperspectral image (HSI) contains hundreds of narrow and contiguous spectral bands, with wavelengths spanning the visible to the infrared spectrum (Chang, 2007). Because of its rich spectral and spatial information, HSI has been widely used in various fields, such as target detection, precision agriculture, mineral exploration, and image fusion (Carrino et al., 2018; Imani & Ghassemian, 2020; Shen et al., 2022; Wei et al., 2019). HSI classification is one of the most vibrant tasks in remote sensing applications and has undoubtedly received a lot of attention, which aims at assigning a spectrum to one certain class (Bioucas-Dias et al., 2013; Camps-Valls et al., 2013; Liu et al., 2020). In the early stages of the study, most methods used exclusively spectral characteristics of a pixel to complete the classification task. Some of these methods focused on feature extraction or dimensionality reduction, such as principal component analysis (PCA; Prasad & Bruce, 2008), independent component analysis (Tu, 2000), and linear discriminant analysis (Bandos et al., 2009). Some other pixel-wise approaches were designed as classifiers, e.g. support vector machines (SVM; Melgani & Bruzzone, 2004), multinomial logistic regression (Pal, 2012), and random forest (Ham et al., 2005).

The problem of purely using spectral information is not only reflected in the lack of spatial information, but also that spectral data is affected by incident illumination, atmospheric effects and instrument noise. These factors greatly limit the accuracy of pixel-level classification methods (Li, Song et al., 2019). To compensate for the inadequacy of single spectral information, some researchers proposed new methods to obtain spatial features for information supplementation (Benediktsson et al., 2005; Camps-Valls et al., 2006; Fauvel et al., 2013; C. Zhao et al., 2017). The results demonstrated the validity of spatial-spectral classification approaches. In addition, some approaches attempted to explore effective feature extraction algorithms for better revealing the structural information of HSI, which can help with subsequent classification tasks. Manifold learning is one of the widely used techniques which can obtain the nonlinear characteristics of data, by combining it with other approaches to improve performance, such as sparse representation (Y. Y. Tang et al., 2014), hypergraph framework (Duan et al., 2021), and geodesic distance (Duan et al., 2022).

In the past few years, with the continuous improvement of computer performance, deep learning has rapidly become a hotspot of research (L. Zhang et al., 2016). The learning process of deep learning methods is fully automated through data-driven technology (Liu, Wu et al. 2022). Due to its outstanding performance, deep learning has been used in object detection, action recognition, and natural language processing (Minar & Naher, 2018). Furthermore, deep learning has also been adopted by many researchers to implement HSI classification tasks. Compared with traditional approaches, deep learning can learn complex and abstract information from data by simulating the recognition process of the brain. Among different deep network models, convolutional neural networks (CNNs) are widely used in HSI classification research(Makantasis et al., 2015, July; Yang et al., 2018). CNNs can capture contextual spatial information in an end-to-end and hierarchical manner. Additionally, CNNs adopt the weight-share mechanism, which greatly reduces the network parameters. Although the role of CNNs in HSI classification is undeniable, it is not good at capturing long-range dependencies of the data.

The advent of the visual Transformer (Dosovitskiy et al., 2020) has given researchers new inspiration to reconsider the image classification process in terms of sequential data to obtain long-range dependencies. The Transformer was first developed by Vaswani et al. (2017) and applied in natural language processing. Numerous experiments have verified the powerful image processing capability of Transformers in the computer vision domain (Carion et al., 2020; Z. Liu et al., 2021; Touvron et al., 2021). In contrast to the Transformer, Tolstikhin et al. (2021) discarded the self-attention operation and proposed an MLP-Mixer (a.k.a., Mixer) architecture based only on multilayer perceptrons (MLPs). It obtains the global perceptual field by matrix transposition and token-mixing projection. Subsequently, Hou et al. (2022) proposed a new MLP-based method named Vision Permutator (ViP), which further splits the information encoding of spatial dimensions to produce more precise positional information. Therefore, ViP has three branches for encoding width, height and channel dimensions, respectively, and it directly sums the features of multiple branches at the end. However, for HSI, the channel information also contains rich features, and aggregating the three branches in one step may result in insufficient fusion of spatial and spectral information.

To address the problems mentioned above, a novel spatial-spectral hierarchical Vision Permutator (S²HViP) network is proposed in this paper for HSI classification. First, in response to the characteristic that HSI feature maps are usually represented in the format of 3-D cubes, we design a hierarchical Vision Permutator (HViP). By utilizing a hierarchical fusion strategy, HViP can not only extract long-range dependencies in parallel from each of the three dimensions of HSI, but also improve the fusion results between dimensions. Then, we design a spectral feature extraction module and a spatial feature extraction module to capture different information. In the spectral module, we group the raw data along the spectral dimension, and the pixels in each group are different spectral tokens. The spectral features are implemented by computing intra- and inter-group spectral correlations. In the spatial module, morphological methods are introduced to facilitate the subsequent network to model deep spatial features. Finally, the features extracted from the different modules are fused to further improve the representation capacity of deep features.

The key contributions of this article can be summarized as follows:

1. We propose a novel S²HViP network for HSI classification. Compared with the mainstream model based on CNNs, S²HViP is good at capturing the global characteristics of data. We design two different modules to extract spectral and spatial long-range dependencies, and these two kinds of global information are aggregated to obtain robust spatial-spectral features for classification.

2. We design an improved ViP for HSI, named HViP. Specifically, it can capture features in the width, height and channel dimensions of HSI, respectively. Based on this, HViP effectively aggregates features in different dimensions through a hierarchical fusion strategy to obtain discrimination features.

Related works

In recent years, HSI classification research has continued to advance. During this period, the researchers have tried a variety of methods to improve classification performance. We divide those methods into traditional approaches, CNN-based approaches, Transformer-based approaches, and MLP-based approaches. In this section, we will review the development of classification methods.

Traditional approaches

In the early time, the spectral characteristics of HSI were the focus of most methods. However, high intra-class spectral variability and low inter-class spectral variability both bring challenges for pixel-wise approaches (Ghamisi et al., 2017; L. He et al., 2018). Considering that the 3-D structure presented by HSI contains a wealth of spectral and spatial features, the introduction of spatial information can effectively improve classification performance (Ghamisi et al., 2018). Morphological profiles (MPs) are generated by morphological transformation approaches, which are used to model spatial information (Pesaresi & Benediktsson, 2001). Based on this, the extended MP is further proposed to finish the HSI classification task (Benediktsson et al., 2005). Random fields and probabilistic graphical models can incorporate spatial features into the classification stage. Markov random field is one of the classic methods used for HSI classification (J. Li et al., 2012), it can be utilized in combination with other technologies, such as SVM (Ghamisi et al., 2014), multinomial logistic regression (Khodadadzadeh et al., 2014), ensemble classifiers (Xia et al., 2015) and active learning (S. Sun et al., 2015). Besides, conditional random field (J. Zhao et al., 2018), superpixel segmentation algorithm (Tu et al., 2020), collaborative representation (Liu et al., 2016) and subspace projection (Wang et al., 2016) are also effective ways to utilize spatial features.

CNN-based approaches

Benefiting from the advancement in parallel computing, deep learning has become the popular method to handle data. Compared with traditional approaches, deep learning-based classifiers use multiple nonlinear networks to hierarchically construct high-level features in an automated manner. Numerous experiments demonstrate its effectiveness in the field of HSI classification (Audebert et al., 2019; Zhu et al., 2017). Typical methods include stacked autoencoder (Chen et al., 2014), deep belief network (T. Li et al., 2014), recurrent neural network (Mou et al., 2017), and CNN (Chen et al., 2016). However, most of the methods in the above models utilize vector inputs except for the CNN, which causes the spatial contextual relationships between pixels to be ignored. The excellent ability of CNN to extract spatial information makes it popular in the HSI field (Guo et al., 2020; X. Li et al., 2019). Paoletti et al. (2019) presented a deep residual network by using pyramidal bottleneck residual blocks. Song et al. (2018) introduced a deep feature fusion network (DFFN) to fuse the outputs of different hierarchical layers to obtain more robust features. In addition to the 2D-CNN-based approaches described above, using a 3-D kernel to extract spectral space features is also a natural solution for HSI classification (Hamida et al., 2018; Zhong et al., 2018). Roy et al. (2020) combined the advantages of 2D-CNN and 3D-CNN to construct a hybrid spectral CNN (HybridSN). The proposed attention mechanism has given a further boost to network performance. Ma et al. (2019) utilized two kinds of branches to form a double-branch multi-attention mechanism network (DBMA) to capture spectral and spatial features, different types of attention mechanisms are used in different branches. Roy et al. (2021) designed a ResNet with spectral attention named attention-based adaptive spectral-spatial kernel ResNet (A²S²K-ResNet). Qing and Liu (2021) presented a multi-scale residual network model with an attention mechanism and extracted spatial-spectral features from different scales. For extracting global spatial information, Zheng et al. (2020) proposed a novel patch-free global learning architecture, which uses the entire HSI as input to a fully convolutional network. Wu et al. (2021) adopted a similar learning architecture and introduced a composite kernel learning network, which extracts spatial-spectral generalized kernel features for accurate classification. Subsequently, a dual-channel convolution network was developed by Yu et al. (2022), it aims to make full use of the global and multi-scale information of HSI.

Transformer-based approaches

ViT (Dosovitskiy et al., 2020) is the first work to construct a vision backbone based exclusively on the Transformer. It demonstrates excellent global feature extraction capabilities. Due to its powerful performance, Transformer is not unexpectedly migrated to HSI classification tasks. J. He et al. (2020) proved that the language model can be applied in HSI classification and proposed HSI-BERT, which is a model consisting of numerous attention layers. Inspired by the HSI-BERT, ZZhong et al. (2022) designed a spectral-spatial Transformer network (SSTN) that replaced convolution operations with attention modules. SSTN was built from several spatial attention modules and spectral correlation modules. The spatial attention module is responsible for capturing the interactions between pixels at all locations, while the spectral correlation module is concerned with the correlation between a compact set of spectral vectors to all locations. Hong et al. (2022) applied the densely sampled method to group the spectral dimensions, and then, a cross-layer Transformer encoder module was employed to learn advanced features from group-wise adjacent bands. Liu, Yu et al. (2022) successfully replaced the traditional convolutional layer with the Transformer, designing a spectral-spatial HSI classification model named DSS-TRM. L. Sun et al. (2022) developed a spectral-spatial feature tokenization Transformer method, which organically combines CNN and Transformer to capture spectral-spatial information and high-level semantic information.

MLP-based approaches

Recently, Tolstikhin et al. (2021) demonstrated that the self-attention layer of Transformers is not necessary and presented a concise alternative model, MLP-Mixer. It is made up of two components: the token-mixing MLP and the channel-mixing MLP. The token-mixing MLP projects feature maps along the spatial dimension to obtain spatial features between different locations. The channel-mixing MLP acts independently on each channel of the feature map to capture the communication between different channels. Inspired by MLP-Mixer, Hou et al. (2022) proposed an effective MLP architecture, Vision Permutator (ViP). The main difference from the MLP-Mixer is that ViP splits spatial feature representation into height-coding and width-coding features respectively, and performs linear projections independently. It allows ViP to capture long-range dependencies along one spatial direction while retaining location information along the other direction. The MLP-based method has also received attention in the hyperspectral field. He and Chen (2021) proposed a pure MLP architecture for the HSI classification task, which demonstrates that MLP networks can provide promising classification performance. Meng et al. (2021) used Mixer as the backbone of the network to extract spatial and spectral information alternately by matrix transposition and MLPs, which allows for interaction between different information. A model based on MLP network and residual learning was presented by X. J. Tang et al. (2022), MLP effectively removes translation invariance and local connectivity. Gong et al. (2022) proposed an MLP architecture, which applies a ladder-like connected structure to obtain contextual interaction for HSI classification tasks. Moreover, Lin et al. (2022) proposed a multi-scale U-shape MLP, which consists of the designed multi-scale channel block and the U-shape MLP.

Proposed framework

In this section, we will introduce the overall structure of S²HViP in detail. Our method contains two modules, the spatial module and the spectral module. Figure 1 shows the framework diagram of our method.

Figure 1. The overall architecture of the proposed S²HViP.

Spectral module

An initial HSI data $I\in {\mathbb{R}}^{H\times W\times B}$ is given, where $H$, $W$, and $B$ represent the height, the width, and the number of spectral bands, respectively. We divide $I$ into non-overlapping $M$ groups along the spectral dimension by 3D-CNN to get grouped data $X_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}}\in {\mathbb{R}}^{H\times W\times M}$, where $M$ is the number of groups, and consider each pixel in the groups as a spectral token. If $B$ is not divisible by $M$, we use the first few bands of $I$ to pad. Then, $X_{\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{u}\mathrm{p}}$ is projected to spectral tokens data $X^{spe}\in {\mathbb{R}}^{H\times W\times C}$ by a linear embedding layer, where $C$ is the hidden layer dimension number. The above process is spectral patch embedding, details of the implementation can be found in Table 1. We enter $X^{spe}$ into the subsequent MLP-based layers. The flowchart of the hierarchical Permutator (HP) block is shown in Figure 2, the biggest difference between our method and Transformer is that we abandon the self-attention operation. A basic HP block consists of LayerNorms, skip connections, hierarchical Permute-MLP (a.k.a., Permute-MLP) and Channel-MLP. We choose the Gaussian error linear unit (GELU; Hendrycks & Gimpel, 2016) as the activation function. The data handling process of HP can be expressed as:

(1) ${\mathit{Y}}^{spe}=\mathrm{P}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{u}\mathrm{t}\mathrm{e}-\mathrm{M}\mathrm{L}\mathrm{P}\left(\mathrm{L}\mathrm{N}\left({\mathit{X}}^{spe}\right)\right)+{\mathit{X}}^{spe}$(1)

(2) ${\mathit{Z}}^{spe}=\mathrm{C}\mathrm{h}\mathrm{a}\mathrm{n}\mathrm{n}\mathrm{e}\mathrm{l}-\mathrm{M}\mathrm{L}\mathrm{P}\left(\mathrm{L}\mathrm{N}\left({\mathit{Y}}^{spe}\right)\right)+{\mathit{Y}}^{spe}$(2)

Figure 2. Flowchart of hierarchical Permutator.

Table 1. Configuration of the spectral module for the Indian Pines dataset.

Layer Name	Input Shape	Kernel Size	Padding	Stride	Filters	Output Shape
Band Padding	(9, 9, 200)	–	–	–	–	(9, 9, 203)
Patch Embedding_1	(1, 9, 9, 203)	(1, 1, 7)	(0, 0, 0)	(1, 1, 7)	1	(1, 9, 9, 29)
Patch Embedding_2	(9, 9, 29)	–	–	–	–	(9, 9, 180)
Hierarchical Permutator × 4	(9, 9, 180)	–	–	–	–	(9, 9, 180)
Global Pooling	(9, 9, 180)	–	–	–	–	(1, 1, 180)

where LN means LayerNorm, ${\mathit{Y}}^{spe}$ and ${\mathit{Z}}^{spe}$ represent the spectral features extracted at different stages respectively. The ${\mathit{Z}}^{spe}$ serves as the input to the subsequent HP blocks.

The visual illustration of the hierarchical Permute-MLP is shown in Figure 3, which has three branches. Linear projections are used to model the input 3-D features separately along their respective dimensions. In the spectral module, we indicate height features, width features and channel features of the extracted feature maps by ${\mathit{X}}_H^{spe}$, ${\mathit{X}}_W^{spe}$, and ${\mathit{X}}_C^{spe}$ respectively. The height and width dimensions represent the spectral correlation within the same group, and the channel dimension stands for the spectral correlation between groups. Specifically, we illustrate the information extraction process using the example of height. We conduct a height-channel permutation operation on ${\mathit{X}}^{spe}$ in the beginning. First, the channel dimension is divided into T parts to yield ${\mathit{X}}_{H_i}^{spe}\in {\mathbb{R}}^{H\times W\times N},i\in \left\{1,\cdots ,T\right\}$, satisfying $C=N\times T$, $N$ is the number of spectra in each part, and it is numerically equal to $H$. Next, the first and third dimensions of ${\mathit{X}}_{H_i}^{spe}$ are permutated to get ${\left({\mathit{X}}_{H_i}^{spe}\right)}^{\mathrm{\top }}\in {\mathbb{R}}^{N\times W\times H}$. All ${\mathit{X}}_{H_i}^{\mathrm{\top }}$ are spliced together along the third dimension in the next step to output ${\left({\mathit{X}}_{H_i}^{spe}\right)}^{\mathrm{\top }}\in {\mathbb{R}}^{N\times W\times C}$. Responsible for encoding ${\left({\mathit{X}}_H^{spe}\right)}^{\mathrm{\top }}$ is a fully connected (FC) layer with weight ${\mathit{W}}_H^{spe}\in {\mathbb{R}}^{C\times C}$. After feature extraction, the original input dimension can be obtained only by performing the inverse process of dimension operation described earlier to yield ${\mathit{X}}_H^{spe}$. In the second branch, we perform a similar operation as above to achieve ${\mathit{X}}_W^{spe}$. The feature encoding of the third branch via a FC layer with weights ${\mathit{W}}_C^{spe}\in {\mathbb{R}}^{C\times C}$ to yield ${\mathit{X}}_C^{spe}$. So far, three types of features $\left({\mathit{X}}_H^{spe},{\mathit{X}}_W^{spe},{\mathit{X}}_C^{spe}\right)$ have been produced by three branches. We use a hierarchical fusion approach to get more recognisable features. The intra-group spectral information is first fused to refine the correlated features between different pixels in the same group. The fused features are then combined with inter-group correlations to obtain spectral long-range dependencies, and inter-group correlations are obtained by mining the spectral features of different groups. To distinguish the importance of different branches, we adopt the split-attention proposed by H. Zhang et al. (2020) to assign weights in the fusion process. The results of the recalibrated fusion of different branches are added to the next FC layer. The fusion operation for correlations of the spectrum can be calculated as follows:

(3) ${\widehat{X}}^{spe}=\text{FC}\left(\text{FC}\left(X_H^{spe}+X_W^{spe}\right)+X_C^{spe}\right)$(3)

Figure 3. Basic structure of the hierarchical Permute-MLP layers.

where $\mathrm{F}\mathrm{C}\left(\cdot \right)$ stands for an FC layer. ${\widehat{X}}^{spe}$ means the output spectral features. After several HP blocks, the spectral information will be sent to the global average pooling layer to obtain the spatial vector for fusion and data classification.

Spatial module

In the spatial module, the extended multi-attribute profile (EMAP) (Dalla Mura, Atli Benediktsson et al., 2010) is utilized to extract preliminary spatial features. Specifically, attribute profiles (APs) (Dalla Mura, Benediktsson et al., 2010) are a multi-level decomposition of the input image based on attribute filters. In contrast to MPs, APs can process images according to a variety of flexibly defined attributes. An AP is obtained by:

(4) $\mathrm{A}\mathrm{P}\left(f\right)=\left\{{\varphi }_n\left(f\right),{\varphi }_{n-1}\left(f\right),\dots ,{\varphi }_1\left(f\right),f,\text{break}{\gamma }_1\left(f\right),\dots ,{\gamma }_{n-1}\left(f\right),{\gamma }_n\left(f\right)\right\}$(4)

where $f$ is the input image, ${\varphi }_j\left(f\right)$ and ${\gamma }_j\left(f\right)$ represent opening and closing operators respectively. Opening and closing are a pair of opposite morphological operators that process images through sliding windows called the structuring element (SE). The size of SE will affect the degree of image processing. Then, we calculate the AP for each principal component (PC) after dimension reduction to get extended APs (EAPs), which can model the needed spatial features.

(5) $\mathrm{E}\mathrm{A}\mathrm{P}=\left\{\mathrm{A}\mathrm{P}\left(\mathrm{P}{\mathrm{C}}_1\right),\mathrm{A}\mathrm{P}\left(\mathrm{P}{\mathrm{C}}_2\right),\dots ,\mathrm{A}\mathrm{P}\left(\mathrm{P}{\mathrm{C}}_{\mathrm{m}}\right)\right\}$(5)

where $\mathrm{m}$ is the number of spectra after PCA. An EMAP contains different EAPs with the aim of extracting different information from the scene, which can be expressed as:

(6) $\mathrm{E}\mathrm{M}\mathrm{A}\mathrm{P}=\left\{\mathrm{E}\mathrm{A}{\mathrm{P}}_{a_1},\mathrm{E}\mathrm{A}{\mathrm{P}}_{a_2},\dots ,\mathrm{E}\mathrm{A}{\mathrm{P}}_{a_k}\right\}$(6)

where $\left\{a_1,a_2,\dots ,a_k\right\}$ denotes the $k$ different attributes.

Mathematically, the input to the spatial module is also $I\in {\mathbb{R}}^{H\times W\times B}$, PCA is first used to process the data to obtain the primary spatial information of the $\mathrm{m}$ PCs. After procedures of the EMAP, we have preliminary spatial features ${\mathit{I}}_{emap}\in {\mathbb{R}}^{H\times W\times G}$, $G$ is the number of the spectral channels. The next step is to divide the obtained features into several tokens evenly and convert the channel dimension number to the hidden layer dimension number. In this module, we choose 2D-CNN to accomplish the spatial token embedding operation. The convolutional layer is composed of $C$ 2-D kernels, The spatial size and stride of each 2-D kernel are both $p$. Then, $C$-dim spatial tokens data ${\mathit{X}}^{spa}\in {\mathbb{R}}^{h\times w\times C}$ are generated, where $h$ and $w$ enumerate the height and the width of tokens, satisfying $h=H/p$ and $w=W/p$. The implementation details of the embedding operation are available in Table 2. These spatial tokens will be sent into a series of HP blocks to extract deeper robust spatial information. Instead of flattening the spatial dimensions, we input the 3-D data ${\mathit{X}}^{spa}$ directly and encode the width and height of the tokens separately. This has the advantage of enabling HViP to capture long-range dependencies in one spatial dimension while preserving position information in the other dimension. In the layer of Permute-MLP, we first perform a weighted fusion of the location information of the height ${\mathit{X}}_H^{spa}$ and width ${\mathit{X}}_W^{spa}$ branches to produce spatial feature representations. Based on the strategy of hierarchical fusion, after fusing the features of ${\mathit{X}}_H^{spa}$ and ${\mathit{X}}_W^{spa}$, the information from the channel branches ${\mathit{X}}_C^{spa}$ is utilized as a complement to further combine to generate the global spatial features. The corresponding calculation process is given as:

(7) ${\widehat{X}}^{spa}=\text{FC}\left(\text{FC}\left(X_H^{spa}+X_W^{spa}\right)+X_C^{spa}\right)$(7)

Table 2. Configuration of the spatial module for the Indian Pines dataset.

Layer Name	Input Shape	Kernel Size	Padding	Stride	Filters	Output Shape
PCA	(9, 9, 200)	–	–	–	–	(9, 9, 3)
EMAP	(9, 9, 3)	–	–	–	–	(9, 9, 114)
Patch Embedding	(9, 9, 114)	(3, 3)	(0, 0)	(3, 3)	180	(3, 3, 180)
Hierarchical Permutator × 4	(3, 3, 180)	–	–	–	–	(3, 3, 180)
Global Pooling	(3, 3, 180)	–	–	–	–	(1, 1, 180)

where ${\widehat{X}}^{spa}$ means the output spatial features. The deep spatial long-range dependencies are obtained after several HP blocks. At the end of the spatial module, the final spatial vector is got by passing through a global average pooling layer, which is also used for fusion and data classification.

Weight fusion

Considering the extracted spectral and spatial features are generated from separate modules, we utilize the weighted fusion approach to obtain spectral-spatial joint features. Specifically, different weighted scores are assigned to spectral features ${\mathit{F}}^{spe}$ and spatial features ${\mathit{F}}^{spa}$, fusing these two features by summing. The fusion progress can be formulated as:

(8) $\mathit{F}=\mathit{\lambda }\times {\mathit{F}}^{spe}+\left(1-\mathit{\lambda }\right)\times {\mathit{F}}^{spa}$(8)

where $F$ denotes the fusion features, $\mathit{\lambda }$ stands the weighting parameter in the range of $\left[0,1\right]$. We send $F$ into a FC layer to get the final feature vector. Note that the length of the vector is the same as the number of classes, and the vector is treated as a class-specific response.

Experiments

Hyperspectral datasets

In this paper, the effectiveness of the proposed method was validated on three benchmark hyperspectral datasets. The detailed data description is presented below.

1. The Indian Pines (IP) was collected by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) sensor over the Indian Pines test site in Northwestern Indiana, USA. It contains 145 × 145 pixels with a spatial resolution of 20 m/pixel (mpp) and 224 spectral bands in the wavelength range from 400 to 2500 nm. We utilize 200 bands after removing four bands containing zero values and 20 noisy bands. The ground truth includes 16 land-cover classes and 10,249 labelled pixels. Figure 4 shows the false-colour composite image and corresponding ground reference of IP. Table 3 outlines the number of labelled samples and classes of IP.

   Figure 4. False-colour composite image and ground-truth map of Indian Pines.

   

   Table 3. Land cover classes and the numbers of samples in the Indian Pines dataset.

   No.	Class	Samples
   1	Alfalfa	46
   2	Corn-notill	1428
   3	Corn-mintill	830
   4	Corn	237
   5	Grass-pasture	483
   6	Grass-trees	7302
   7	Grass-pasture-mowed	28
   8	Hay-windrowed	478
   9	Oats	20
   10	Soybean-notill	972
   11	Soybean-mintill	2455
   12	Soybean-clean	593
   13	Wheat	205
   14	Woods	1265
   15	Buildings-Grass-Trees-Drives	386
   16	Stone-Steel-Towers	93

• (2) The Salinas Valley (SV) was recorded by the AVIRIS sensor over Salinas Valley, California, USA. The available dataset is composed of 204 bands with 512 × 217 pixels with a spatial resolution of 3.7 m after the low signal-to-noise ratio (SNR) bands were removed. The ground truth includes 16 land-cover classes and 54,129 labelled pixels. Figure 5 shows the false-colour composite image and corresponding ground reference of SV and the details of labelled samples and classes are presented in Table 4.

  Figure 5. False-colour composite image and ground-truth map of Salinas Valley.

  

  Table 4. Land cover classes and the numbers of samples in the Salinas Valley dataset.

  No.	Class	Samples
  1	Broccoli-green-weeds-1	2009
  2	Broccoli-green-weeds-2	3726
  3	Fallow	1976
  4	Fallow-rough-plow	1394
  5	Fallow-smooth	2678
  6	Stubble	3959
  7	Celery	3579
  8	Grapes-untrained	11,271
  9	Soil-vinyard-develop	6203
  10	Corn-senesced-green-weeds	3278
  11	Lettuce-romaine-4 wk	1068
  12	Lettuce-romaine-5 wk	1927
  13	Lettuce-romaine-6 wk	916
  14	Lettuce-romaine-7 wk	1070
  15	Vinyard-untrained	7268
  16	Vinyard-vertical-trellis	1807

• (3) The University of Pavia (UP) was acquired by the ROSIS sensor during a flight campaign over Pavia, Northern Italy. It consists of 610 × 340 pixels with a spatial resolution of 1.3 m and has 115 spectral channels in the wavelength range from 0.43 to 0.86 µm. We utilize 103 bands after removing 12 noisy bands. The ground truth includes 9 land-cover classes and 42,776 labelled pixels. Figure 6 shows the false-colour composite image and corresponding ground reference of UP. The detailed numbers of labelled samples and classes of UP are listed in Table 5.

  Figure 6. False-colour composite image and ground-truth map of University of Pavia.

  

  Table 5. Land cover classes and the numbers of samples in the University of Pavia dataset.

  No.	Class	Samples
  1	Asphalt	6631
  2	Meadow	18,649
  3	Gravel	2099
  4	Trees	3064
  5	Metal sheets	1345
  6	Bare soil	5029
  7	Bitumen	1330
  8	Bricks	3682
  9	Shadows	947

Experimental configuration

To verify the effectiveness of the proposed method, we compare it with several recently proposed HSI classification methods, including DFFN (song et al., 2018), SSRN (Zhong et al., 2018), DBMA (Ma et al., 2019), pResNet (Paoletti et al., 2019), HybridSN (Roy et al., 2020), A²S²K-ResNet (Roy et al., 2021), FreeNet (Zheng et al., 2020), DPSCN (Dang et al., 2021), SSTN (ZZhong et al., 2022), and Mixer (Tolstikhin et al., 2021). For the fairness of the experiment, we adopt unified measurements to illustrate the effectiveness of our framework. These measurements are overall accuracy (OA), average accuracy (AA) and Kappa coefficient (Kappa).

Details of experimental configurations are presented as follows. For the IP dataset, 1% of labelled samples are randomly chosen from each class to form the training set, and the rest are used as the testing set. For UP and SV datasets, 0.5% and 99.5% of labelled samples are selected for the training set and the testing set. The stochastic gradient descent (SGD) is utilized to train parameters in the whole architecture. The batch size is 50 and the learning rate is 0.001. We set the size of the input patch to 9. And the network runs 200 epochs at a time to get the final result. In the spatial module, two kinds of APs are selected to calculate EMAPs: the moment of inertia and area. For the moment of inertia, all three datasets are computed with the same values {0.2, 0.3, 0.4, 0.5}. The difference is the attribute of area, which has a total of 14 values, different datasets have different values. For the IP dataset, the initial value of the area value is set to 100 with a step size of 400. For the SV dataset, the values of the area are from the range of 270 to 7300 in the step of 540. For the UP dataset, the area values ranged from 770 to 10,769 and the step size is being 769.

Results analysis

For the aforementioned classification approaches, Tables 6–8 show the OAs, AAs, and Kappas on IP, SV, and UP datasets, respectively. As can be seen, the proposed method performs better than other methods in three datasets. In the following content, we will further analyse the experimental results of three datasets.

Table 6. Classification results for the Indian Pines dataset using 1% training samples.

No.	DFFN	SSRN	DBMA	pResNet	Hybrid-SN	FreeNet	A^2S^2K-ResNet	DPSCN	SSTN	Mixer	Proposed
1	79.77	94.65	95.58	84.19	75.11	84.04	99.07	100.0	95.12	93.72	93.26
2	61.77	79.02	68.42	57.79	59.63	71.54	69.62	74.55	80.47	63.29	68.04
3	40.40	66.74	61.91	33.19	46.20	63.81	61.27	64.66	67.01	57.98	71.85
4	36.28	0.06	74.70	24.23	34.70	85.84	48.50	53.89	48.68	48.25	75.26
5	47.03	70.15	69.08	57.35	67.67	87.95	77.47	81.98	68.91	68.77	77.20
6	83.80	98.99	89.85	84.92	92.23	92.58	94.61	98.73	98.53	92.31	95.10
7	94.40	96.00	100.0	85.20	88.00	48.40	100.0	100.0	100.0	90.00	100.0
8	89.03	98.63	95.00	89.58	74.03	100.0	95.49	92.09	99.83	99.47	100.0
9	97.05	52.94	94.12	94.71	87.65	28.82	100.0	100.0	96.47	98.24	100.0
10	68.00	74.61	77.47	67.91	57.81	62.48	74.74	76.31	72.99	76.48	78.96
11	69.74	88.22	81.72	66.54	73.53	74.33	83.09	84.75	86.41	91.65	93.93
12	33.35	31.94	63.20	27.36	39.64	74.23	49.59	69.68	75.88	63.08	88.30
13	86.04	86.29	88.42	82.72	61.63	99.43	80.05	90.99	87.43	86.44	92.22
14	82.69	99.78	96.73	85.61	89.10	84.30	99.03	99.11	98.75	99.77	100.0
15	44.70	22.64	50.73	42.95	29.24	89.30	57.44	59.79	53.13	83.71	91.04
16	95.44	99.67	94.00	96.00	56.33	94.12	99.33	99.89	98.22	99.78	96.44
OA	65.39	78.15	78.27	63.67	66.03	76.94	78.27	81.74	82.40	80.89	86.58
AA	69.34	72.91	81.31	67.52	64.53	77.58	80.58	84.16	82.99	82.06	88.85
Kappa	60.68	74.76	75.23	58.57	61.07	73.51	74.98	79.10	79.82	78.14	84.70
Train(s)	14.86	138.98	67.00	16.04	10.73	4.15	140.55	6.63	11.73	17.15	27.39
Test (s)	2.91	2.98	13.30	3.51	4.07	0.01	3.57	1.17	1.76	3.29	5.04

Table 7. Classification results for the Salinas Valley dataset using 0.5% training samples.

No.	DFFN	SSRN	DBMA	pResNet	Hybrid-SN	FreeNet	A^2S^2K-ResNet	DPSCN	SSTN	Mixer	Proposed
1	90.79	99.83	99.88	87.30	99.23	100.0	98.92	98.81	98.53	100.0	100.0
2	98.81	100.0	99.96	99.98	99.97	90.90	99.84	100.0	100.0	100.0	100.0
3	76.41	99.17	95.49	93.25	98.37	99.87	99.37	99.98	93.50	97.85	99.92
4	99.15	99.82	99.04	99.67	99.20	99.04	99.97	89.94	99.56	99.96	99.98
5	96.02	98.59	95.67	97.77	96.96	99.05	96.84	97.63	95.74	93.95	97.04
6	98.29	99.41	98.71	98.12	98.21	99.42	98.58	99.13	99.12	99.39	98.95
7	99.81	100.0	99.99	99.71	98.48	99.03	100.0	100.0	99.98	100.0	100.0
8	76.49	88.15	89.39	77.66	92.79	99.48	83.93	84.32	77.59	95.87	97.14
9	99.26	100.0	100.0	99.62	99.94	99.01	100.0	99.99	100.0	100.0	100.0
10	93.67	98.22	97.46	93.73	98.14	95.19	98.45	97.90	96.91	98.98	99.45
11	81.66	98.79	81.72	96.75	94.90	94.59	99.72	89.59	98.74	99.70	99.88
12	99.00	100.0	99.57	97.62	97.64	99.94	99.95	99.98	99.44	99.72	99.99
13	99.33	87.49	98.37	99.08	98.90	100.0	98.29	99.16	98.08	98.22	99.76
14	98.51	98.54	97.65	98.98	91.02	99.98	98.98	98.70	98.55	99.47	99.00
15	85.91	87.47	87.80	69.81	93.65	96.28	89.03	93.03	90.33	95.27	97.99
16	96.01	99.77	99.08	97.66	99.48	99.81	99.58	99.45	99.39	100.0	100.0
OA	90.53	95.31	95.68	89.61	96.64	97.98	94.68	94.92	93.14	97.96	98.85
AA	93.07	97.20	97.32	94.17	97.31	98.23	97.59	96.73	96.59	98.65	99.32
Kappa	89.50	94.79	94.83	88.46	96.26	97.75	94.09	94.36	92.39	97.73	98.72
Train(s)	28.21	335.34	106.25	25.83	11.67	21.06	340.16	15.49	23.10	34.59	62.32
Test (s)	13.87	16.46	43.76	14.18	10.63	0.06	19.65	6.50	9.88	17.64	28.85

Table 8. Classification results for the University of Pavia dataset using 0.5% training samples.

No.	DFFN	SSRN	DBMA	pResNet	Hybrid-SN	FreeNet	A^2S^2K-ResNet	DPSCN	SSTN	Mixer	Proposed
1	92.65	98.42	97.24	88.78	81.54	94.32	96.62	98.50	98.35	99.96	100.0
2	97.24	99.92	99.65	97.58	99.64	96.94	99.79	99.93	99.80	99.97	99.92
3	47.24	83.12	82.92	60.78	66.42	81.18	84.22	97.26	89.80	89.78	90.61
4	87.21	87.84	89.16	87.70	61.33	99.09	91.03	86.92	84.73	82.22	87.94
5	98.37	99.99	99.57	99.19	72.02	99.89	99.81	99.96	99.87	99.97	100.0
6	74.56	94.34	94.92	65.37	86.11	92.23	87.75	93.65	96.39	99.00	99.96
7	64.09	80.75	81.08	59.82	74.26	87.16	89.06	95.22	95.38	99.38	99.86
8	85.42	92.85	94.61	73.03	63.77	88.47	92.77	93.43	87.65	97.48	98.54
9	96.01	98.52	92.23	98.61	40.37	94.02	99.17	99.17	96.46	94.93	95.16
OA	88.65	96.11	95.97	86.70	84.81	94.27	95.05	97.19	96.35	97.74	98.40
AA	82.53	92.86	92.38	81.21	71.72	92.56	90.83	96.01	94.27	95.85	96.89
Kappa	84.76	94.80	94.63	82.10	79.45	92.38	93.38	96.24	95.13	96.98	97.87
Train(s)	28.31	258.12	46.61	25.24	35.59	20.32	274.21	10.82	18.65	28.49	54.39
Test (s)	13.60	7.33	29.05	13.05	25.51	0.10	11.09	4.17	7.25	13.46	38.02

As shown in Table 6, improvements in terms of OA compared to DFFN, SSRN, DBMA, pResNet, HybridSN, FreeNet, A²S²K-ResNet, DPSCN, SSTN, Mixer are 21.19%, 8.43%, 8.31%, 22.91%, 20.55%, 9.64%, 8.31%, 4.84%, 4.18%, and 5.69% respectively. In the case of the IP dataset, the diversity of land cover and the imbalance between sample sizes pose a challenge to the classification of this dataset, especially the quantitative difference between categories 7 and 9 and other categories. Even so, compared with the previous approaches, the proposed S²HViP obtains the best AA and achieves satisfactory accuracy in the class with a small number of samples. In comparison with the second-best method DPSCN, the improved AA value is 4.69%. This is probably due to the great power of capturing long-range spatial interactions of our method, resulting in a lower error score rate when distinguishing between different land-cover classes. The classification map of the IP dataset is displayed in Figure 7. It is noteworthy that our delineated boundaries are sharper and straighter in the boundary regions of different feature classes, which also reflects the fact that S²HViP can better extract global distribution features. This helps to generate a more reasonable regional division.

Figure 7. Classification maps resulting from different methods for Indian Pines dataset. (a) Ground truth. (b) DFFN (OA:65.39%). (c) SSRN (OA:78.15%). (d) DBMA (OA:78.27%). (e) pResNet (OA:63.67%). (f) HybridSN (OA:66.03%). (g) FreeNet (OA:76.94%). (h) A²S²K-ResNet (OA:78.27%). (i) DPSCN (81.74%). (j) SSTN (OA:82.40%). (k) Mixer (OA:80.89%). (l) Proposed (OA:86.58%).

Table 7 shows the classification results of the SV dataset. The characteristic of this dataset is that similar objects are clustered in a similar area. Therefore, the SV dataset has very strong spatial similarity and rich texture features. Since FreeNet uses the whole HSI as input, the spatial features are well preserved. Therefore, FreeNet gets good results of 97.98% OA. However, the proposed S²HViP obtains a better OA of 98.85%. This can be attributed to the long-distance dependencies extracted by S²HViP, which may be important for classification tasks in large-scale regions. Compared with the MLP-based method Mixer, our approach shows 0.66%, 1.04%, and 0.98% improvements in terms of OA, AA, and Kappa. One reason is that S²HViP further splits the spatial information into width and height dimensions to encode the information separately, thus obtaining finer location information. Another reason is the distribution fusion strategy facilitates better interaction of the information in the three dimensions of the feature map. The classification map of the dataset is shown in Figure 8. The map still demonstrates strong classification capabilities of S²HViP. Our method has fewer misclassification points in comparison with other methods, which means our classification graph is closer to the truth of the SV dataset.

Figure 8. Classification maps resulting from different methods for the Salinas Valley dataset. (a) Ground truth. (b) DFFN (OA:90.53%). (c) SSRN (OA:95.31%). (d) DBMA (OA:95.68%). (e) pResNet (OA:89.61%). (f) HybridSN (OA:96.64%). (g) FreeNet (OA:97.98%). (h) A²S²K-ResNet (OA:94.68%). (i) DPSCN (94.92%). (j) SSTN (OA:93.14%). (k) Mixer (OA:97.96%). (l) Proposed (OA:98.85%).

Focusing on the UP dataset, Table 8 shows the classification results. As shown in Table 8, the proposed approach achieves the best OA, AA, and Kappa values of 98.40%, 96.89%, and 97.87%, respectively, while the SSTN network achieves 96.35%, 94.27%, and 95.13%, respectively. The promotions of the OA, AA and Kappa values are 2.05%, 2.62%, and 2.74%, respectively. Compared to Transformer-based methods, S²HViP abandons the self-attention operation. Nevertheless, the results strongly prove the effectiveness of the MLP-based structure. According to Figure 9, compared to the classification maps of other methods, our method has fewer pixels of other categories appearing within the same category region. On the one hand, our method can efficiently extract long-range dependencies through global receptive fields, reducing the likelihood that pixels within connected regions are mistaken as other high-similarity objects. On the other hand, the introduction of the morphological approaches can complement local spatial features, which helps the subsequent information extraction process. As a result, S²HViP obtains competitive results in comparison with other methods.

Figure 9. Classification maps resulting from different methods for the University of Pavia dataset. (a) Ground truth. (b) DFFN (OA:88.65%). (c) SSRN (OA:96.11%). (d) DBMA (OA:95.97%). (e) pResNet (OA:86.70%). (f) HybridSN (OA:84.81%). (g) FreeNet (OA:94.27%). (h) A²S²K-ResNet (OA:95.05%). (i) DPSCN (97.19%). (j) SSTN (OA:96.35%). (k) Mixer (OA:97.74%). (l) Proposed (OA:98.40%).

Investigation of different numbers of training samples

To further verify the effectiveness of S²HViP, we investigate the performance of the different approaches with different training sample sizes. For IP and SV datasets, 0.1%, 0.3%, 0.5%, 1% and 3% of samples are randomly selected. And for the UP dataset, we randomly choose 0.3%, 0.5%, 1%, 3% and 5% of samples as training sets. The specific results counted on different datasets are displayed in Figure 10. Overall, the proposed S²HViP shows better performance. As the sample size increases, the performance of all methods tends to improve in general, and the differences between them become smaller. Nevertheless, when the sample is insufficient, our method still shows better performance than other methods. The results in the figure demonstrate the robustness and efficiency of S²HViP, which could be attributed to its excellent long-range dependencies extraction capability. In the case of limited training samples, our method is still able to capture satisfactory global distribution information, resulting in better classification results.

Figure 10. The comparison results using different training sample ratios on (a) Indian Pines, (b) Salinas Valley, (c) University of Pavia.

Investigation of different sizes of patch

To explore the influence of different patch sizes on the experimental results, on the basis of 9 × 9, we further conducted experiments on the patch sizes of 3 × 3, 6 × 6 and 12 × 12, and the results are shown in Table 9. For IP and UP datasets, OA first increases and then decreases as the patch size increases. In addition, S²HViP makes good use of the spatial aggregation properties of the SV dataset. The results show that the value of OA continues to improve as the patch size increases. This indicates that the SV dataset is more adaptive to large patch sizes. Considering the running time of the model, we finally choose 9 × 9 as the consistent patch size of the experiment.

Table 9. Overall accuracy(OA) on three datasets with different patch sizes.

	3 × 3	6 × 6	9 × 9	12 × 12
IP	75.26	83.15	86.58	85.27
SV	97.08	98.54	98.85	99.03
UP	97.61	98.19	98.40	97.84

Ablation study

The proposed S²HViP consists of a spatial module and a spectral module. The purpose of the former is to obtain robust spatial features. The latter is dedicated to extracting spectral information from the rich spectral features. To investigate the respective functions of these two parts, we divide the network into two individual modules for separate experiments and observe the respective classification results. In other words, we use either the spatial module or the spectral module alone as the overall feature extraction framework for the network. For the spatial module, we further explore the significance of the EMAP module. As can be gathered from Table 10, the spectral module outperforms the spatial module overall without the initial modelling of the image using EMAP. After modelling the spatial information using EMAP, the resulting spatial features are more easily extracted in subsequent modules. The sliding window mechanism of SE also helps to compensate for the local spatial information of feature maps. These factors lead to a significant improvement in the performance of the spatial module. In addition, the global spatial-spectral information obtained by fusing the different features gives the best results when using dual modules.

Table 10. The ablation analysis of different modules (OA%) for three datasets.

	IP	SV	UP
Spectral module	75.32	95.30	94.99
Spatial module w/o EMAP	76.13	94.94	94.59
Spatial module	85.85	98.76	98.05
Dual modules	86.58	98.85	98.40

To further validate the effectiveness of the hierarchical fusion strategy, we compare the performance gap between HViP and direct weighted fusion approaches, i.e. ViP. The results are listed in Table 11. On three different datasets, the HViP achieves 1.78% (IP), 0.17% (SV), and 0.45% (UP) improvement in terms of OA compared to ViP. The experimental results demonstrate the effectiveness of the hierarchical fusion strategy, which can help the three dimensions of the feature map to fuse better and generate more discriminative global features.

Table 11. The ablation analysis of hierarchical integration strategy (OA%) for three datasets.

	IP	SV	UP
ViP	84.80	98.68	97.95
HViP	86.58	98.85	98.40

Conclusion

In this article, we present a novel MLP-based deep classification method. Considering that HSI is presented as 3-D data, we encode the three dimensions individually and combine them through a hierarchical fusion strategy. The strategy is demonstrated to be effective in improving the ability of feature extraction. The proposed S²HViP contains two feature extraction modules for learning the spectral and spatial information of HSI respectively. In the spectral module, we first capture intra- and inter-group spectral correlations by grouping. Then, we aggregate different correlations to refine long-range spectral dependencies. In the spatial module, we introduce morphological approaches to better model spatial features. Based on this, deep spatial information is further captured through MLPs. Finally, the information interaction between the spectral and spatial domains is achieved by weighted fusion, which effectively enhances the classification results. Experiments on three benchmark HSI datasets demonstrate the satisfactory performance of S²HViP and show the potential of the MLP-based network for HSI classification.

In the future, we will consider the application of more innovative morphological methods to the improve feature extraction capability of the network. Moreover, the proposed method still has room for improvement in terms of time consumption. So we will also explore a lightweight MLP-based model for classification.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study are available at http://www.ehu.eus/ccwintco/index.php?%20title=Hyperspectral_Remote_Sensing_Scenes.

Additional information

Funding

This work was supported in part by the National Natural Science Foundation of China under Grant No. 62071204 and by the Natural Science Foundation of Jiangsu Province under Grant No. BK20201338.