# Superpixel Graph Contrastive Clustering with Semantic-Invariant
Augmentations for Hyperspectral Images

Jianhan Qi, Yuheng Jia, Hui Liu, Junhui Hou Jianhan Qi is with the School of
Software Engineering, Southeast University, Nanjing 210096, China (e-mail:
qjh@seu.edu.cn).Yuheng Jia is with the School of Computer Science and
Engineering, Southeast University, Nanjing 210096, China, and also with the
Key Laboratory of New Generation Artificial Intelligence Technology and Its
Interdisciplinary Applications (Southeast University), Ministry of Education,
China (e-mail: yhjia@seu.edu.cn).Hui Liu is with the School of Computing
Information Sciences, Caritas Institute of Higher Education, Hong Kong
(e-mail: hliu99-c@my.cityu.edu.hk).Junhui Hou is with the Department of
Computer Science, City University of Hong Kong, Hong Kong (e-mail:
jh.hou@cityu.edu.hk). 0009-0004-3810-083X 0000-0002-3907-6550
0000-0003-2159-025X 0000-0003-3431-2021

###### Abstract

Hyperspectral images (HSI) clustering is an important but challenging task.
The state-of-the-art (SOTA) methods usually rely on superpixels, however, they
do not fully utilize the spatial and spectral information in HSI 3-D
structure, and their optimization targets are not clustering-oriented. In this
work, we first use 3-D and 2-D hybrid convolutional neural networks to extract
the high-order spatial and spectral features of HSI through pre-training, and
then design a superpixel graph contrastive clustering (SPGCC) model to learn
discriminative superpixel representations. Reasonable augmented views are
crucial for contrastive clustering, and conventional contrastive learning may
hurt the cluster structure since different samples are pushed away in the
embedding space even if they belong to the same class. In SPGCC, we design two
semantic-invariant data augmentations for HSI superpixels: pixel sampling
augmentation and model weight augmentation. Then sample-level alignment and
clustering-center-level contrast are performed for better intra-class
similarity and inter-class dissimilarity of superpixel embeddings. We perform
clustering and network optimization alternatively. Experimental results on
several HSI datasets verify the advantages of the proposed method, e.g., on
India Pines, our model improves the clustering accuracy from 58.79% to 67.59%
compared to the SOTA method.

###### Index Terms:

Hyperspectral image, superpixel, graph learning, contrastive learning,
clustering

##  I Introduction

Hyperspectral remote sensing techniques employ sensors to collect the
reflectance of land-cover materials in hundreds of narrow and contiguous
spectral bands, generating hyperspectral images (HSI), which have been widely
applied in various fields, such as agricultural monitoring[1], mineral
exploration[2] and urban planning[3]. Grouping HSI into several classes is
essential in many applications[4, 5, 6, 7]. However, due to the high cost of
manual annotating, obtaining pixel-level labels is difficult.

Unsupervised clustering methods[8, 9, 10, 11] are gaining increasing
attention. These methods aim to mine relationships between pixels in the
absence of labels and assign similar pixels to the same cluster. For example,
center-based clustering methods, e.g., K-means[12] and fuzzy c-means[13],
assume that instances have a spherical cluster structure, which may be not
satisfied in the raw feature space of HSI. Subspace-based clustering methods,
e.g., sparse subspace clustering[14], low-rank subspace clustering[15], and
dense subspace clustering[16], divide instances into different low-dimensional
subspaces, which are susceptible to redundant spectral features. Additionally,
considering the large quantity of HSI pixels, subspace clustering methods may
lead to unbearable computational costs because they need to calculate the
self-represented coefficient matrix with time complexity of
$\mathcal{O}(n^{3})$, where $n$ is the number of pixels. Therefore, some
competitive subspace-based methods are limited to experimenting on a subset of
the complete HSI[17, 18].

To this end, superpixel segmentation techniques, e.g., ERS[19] and SLIC[20],
which divide HSI into irregular superpixels with high homogeneity, have been
widely used in HSI clustering. In these methods, pixels within the same
superpixel are assumed to belong to the same class, and superpixel-level
clustering is performed. As the number of superpixels is much smaller than
that of pixels, the computational efficiency can be greatly improved by
transforming the clustering task from pixel level to superpixel level.
Moreover, superpixel representations are smoother and less susceptible to
spectral noise[21]. Some methods have combined superpixel segmentation with
traditional clustering techniques, resulting in competitive performance. For
example, Hinojosa et al. took advantage of HSI’s neighboring spatial
information to improve the superpixel-level subspace clustering algorithm[22].
Zhao et al. introduced both global and local similarity matrices to boost the
clustering accuracy[8].

Since superpixels have a clear spatial adjacent relationship with the neighbor
superpixels, it is natural to process superpixels as graph nodes with edges
connected to their spatial neighbors, and many efforts have been devoted to
extracting structured information of superpixels with graph learning methods,
e.g., graph convolutional networks (GCN)[23] and graph attention networks
(GAT)[24]. The most popular training framework, graph autoencoder (GAE)[25],
takes the reconstruction of adjacency matrices as the optimization target and
learns low-dimensional representations of graph nodes. Based on GAE, Zhang et
al. achieved HSI clustering by building spatial and spectral similarity graphs
of superpixels and extracting both spatial and spectral features
simultaneously with dual graph convolutional encoders[10]. Ding et al.
proposed a low-pass GCN with the layer-wise attention mechanism to extract the
smooth and local features of superpixels[11]. Moreover, the auxiliary
distribution loss is often introduced to guide the more concentrated cluster
structure[26, 27]. Generally, such methods model the relationships between
superpixels with a similarity matrix and gain superpixel representations by
GAE. However, they mainly have two drawbacks. First, the input superpixel
features are directly represented as the average of internal pixels, ignoring
the 3-D structure of HSI at the pixel level, thus not fully utilizing the
spatial and spectral information of HSI. Second, the learned superpixel
representations of the above methods may be not clustering-friendly due to
their training objectives are not clustering-oriented. Specifically, features
that are easy to cluster tend to exhibit intra-class similarity and inter-
class dissimilarity, but GAE-based methods mainly aim to reconstruct the input
graph structure, thus limiting the clustering performance.

To overcome the aforementioned drawbacks, we propose a superpixel-level
clustering method for HSI with pixel-level pre-training, named superpixel
graph contrastive clustering (SPGCC). First, we exploit the 3-D structure of
HSI by dividing the original HSI into 3-D cubes. Then, 3-D and 2-D hybrid
convolutions are performed on pixel cubes to extract the high-order spatial
and spectral features of HSI simultaneously which can be preserved in
superpixel-level representations. Next, we aggregate the topological
information of superpixels from neighborhoods by GCN. To learn the
discriminative and clustering-friendly superpixel representations, we utilize
contrastive learning (CL), which brings samples close to their augmented views
in the embedding space while pushing different samples’ representations far
apart[28, 29, 30]. Proper data augmentations are crucial for CL, as augmented
views and original samples should share the same semantics. To this end, we
design two effective data augmentations for superpixels, i.e., pixel sampling
augmentation and model weight augmentation, to obtain semantic-invariant
augmented views. Additionally, traditional CL-based clustering methods treat
different samples as negative pairs and push their representations to be far
apart from each other[31], even for samples in the same class, which is
harmful to clustering accuracy. Differently, we perform sample-level alignment
and clustering-center-level contrast on superpixels to obtain better intra-
class similarity and inter-class dissimilarity. To mitigate the impact of
outliers, we select samples with high confidence to recompute clustering
centers. In the whole framework, we alternatively conduct clustering and
network optimization, using clustering results to guide the optimization. The
proposed SPGCC has an improvement of 8.8%, 2.4%, and 3.2% over the best
compared methods on three different HSI datasets.

The main contributions of this paper are summarized as follows:

  1. 1)

We perform pixel-level pre-training to extract high-order spatial and spectral
information of HSI, where the pre-training network is completely separated
from the subsequent clustering network and can be applied to large-scale HSI.

  2. 2)

We propose a novel superpixel graph contrastive clustering (SPGCC) method and
design two semantic-invariant data augmentations for superpixels, i.e., model
weight augmentation and pixel sampling augmentation, which can obtain reliable
positive samples and improve contrastive learning.

  3. 3)

SPGCC utilizes a more reasonable optimization target compared with traditional
contrastive clustering methods and recomputes high-confidence clustering
centers considering the outliers. Experimental results on three HSI benchmark
datasets demonstrate the effectiveness of our proposed method.

The rest of this paper is organized as follows. Section II briefly reviews
preliminary about HSI pixel convolution and graph contrastive learning.
Section III presents the proposed SPGCC algorithm in detail, including pixel-
level pre-training, graph construction, semantic-invariant superpixel
augmentations, and graph contrastive clustering. In Section IV, we provide
experimental results and analyses on several HSI datasets. Finally, we
summarize our work in Section V.

##  II Preliminary

###  II-A HSI Pixel Convolution

Convolutional neural networks (CNNs) are powerful tools for extracting high-
order features from images. For HSI, a common idea is to expand the pixel into
a square neighborhood centered on itself, forming pixel cubes. Conventional
2-D CNNs simply treat the pixel cube as a multi-channel image and perform
convolution at each channel, ignoring the rich spectral features of HSI.
Therefore, many efforts have been devoted to extracting spatial and spectral
features from HSI simultaneously with 3-D CNNs[32, 33, 34]. Considering the
high computational cost of 3-D CNNs, Roy et al. proposed 3-D and 2-D hybrid
CNNs to accelerate computation efficiency[35].

Given an HSI cube and a 3-D convolutional kernel, the kernel moves along the
width, height, and depth (spectral) dimensions respectively with a fixed step.
At each position, the convolution is defined as the Hadamard product of the
kernel and the related part of the HSI cube. Similar to 2-D convolution, 3-D
convolution is also conducted in a multi-channel way to extract features of
different scales. The output value in the $j$-th channel of the $i$-th layer
at spatial position $(x,y,z)$, denoted as $v_{i,j}^{x,y,z}$, is calculated by:

$\displaystyle v_{i,j}^{x,y,z}=$ (1)
$\displaystyle\phi(\sum\limits_{\tau=1}^{d_{l-1}}\sum\limits_{\lambda=-\eta}^{\eta}\sum\limits_{\rho=-\gamma}^{\gamma}\sum\limits_{\sigma=-\delta}^{\delta}\omega_{i,j,\tau}^{\sigma,\rho,\lambda}\times
v_{i-1,\tau}^{x+\sigma,y+\rho,z+\lambda}+b_{i,j}),$

where $d_{l-1}$ is the number of channels in the $(l-1)$-th layer,
$\omega_{i,j,\tau}\in\mathbb{R}^{(2\delta+1)\times(2\gamma+1)\times({2\eta+1})}$
is the kernel’s weight parameter in the $j$-th channel of the $i$-th layer
related to the $\tau$-th channel of the last layer, superscripts
$(\sigma,\rho,\lambda)$ are indices along corresponding dimensions, $b_{i,j}$
is the bias, and $\phi$ is the activation function.

HSI pixel convolution has been widely used in classification tasks,
significantly improving the accuracy, especially with 3-D convolution. For
clustering tasks, although 2-D convolution has been utilized as the feature
extraction network[9, 36], 3-D convolution is much less used. An intuitive
reason is that 3-D convolution has higher complexity and some clustering
networks require a complete batch input, resulting in excessive GPU memory
consumption.

Figure 1: The diagram of superpixel graph contrastive clustering (SPGCC).
First, the input HSI is divided into several 3-D pixel cubes, and superpixels
are obtained by entropy rate superpixel segmentation (ERS). At pixel-level
pre-training, high-order spatial and spectral features are extracted by 3-D
and 2-D hybrid CNNs encoder and decoder. At superpixel-level contrastive
clustering, model weight augmentation (MWA) and pixel sampling augmentation
(PSA) are combined to generate semantic-invariant augmentations, i.e.,
unshared dual branches in the last layer of graph encoder output different
augmented views for superpixels and sampled internal pixels, respectively.
Then, K-means is performed to get clustering results, and high-confidence
clustering centers of two augmented superpixel views
$\mathbf{c}^{1},\mathbf{c}^{2}$ are recomputed based on samples close to the
original centers. Finally, the network is trained with sample-level alignment
(SLA) and clustering-center-level contrast (CLC), where the different views
are encouraged to be similar at both the sample level and clustering center
level while the contrast is only performed at the clustering center level.
Clustering and network optimization are performed alternatively.

###  II-B Graph Contrastive Learning

Several works have applied contrastive learning for graph data due to its
remarkable representation ability[37, 38, 39, 40, 41]. Here, we focus on node-
level graph contrastive learning (GCL) and divide it into three parts: graph
data augmentation, graph encoder network, and contrastive loss function. A
graph is denoted as $g=(\mathcal{V},\mathcal{E})$ where $\mathcal{V}$ is the
set of $N$ nodes and $\mathcal{E}$ is the set of $M$ edges. The graph
structure is represented by an adjacency matrix
$\mathbf{A}\in\mathbb{R}^{N\times N}$ and node attributes are denoted as
$\mathbf{X}\in\mathbb{R}^{N\times d}$. We use $\Gamma$ to denote a set of
different graph data augmentations and then augmented views are defined as:

$\mathbf{A}^{k},\mathbf{X}^{k}=\gamma^{k}(\mathbf{A},\mathbf{X});\gamma^{k}\in\Gamma.$
(2)

Node embeddings are obtained through graph encoder network $f(\cdot)$, e.g.,
graph convolutional network[23] or graph attention network[24].

Given two different views of graph nodes, their embeddings are calculated as
$\mathbf{Z}^{1}=f(\mathbf{A}^{1},\mathbf{X}^{1}),\mathbf{Z}^{2}=f(\mathbf{A}^{2},\mathbf{X}^{2})$
and $\mathbf{z}_{i}^{k}(k=1,2)$ is the $i$-th node’s embedding from the $k$-th
view. Most GCL methods follow the InfoMax[42] principle, which maximizes the
mutual information between two augmented views of the graph. According to the
most commonly used InfoNCE loss[29], embeddings of the same node from
different views form positive sample pairs, while others form negative sample
pairs. Finally, pairwise objective for
$(\mathbf{z}_{i}^{1},\mathbf{z}_{i}^{2})$ is defined as follows:

$\mathcal{L}(\mathbf{z}_{i}^{1},\mathbf{z}_{i}^{2})=-log\frac{exp(sim(\mathbf{z}_{i}^{1},\mathbf{z}_{i}^{2})/\tau)}{exp(sim(\mathbf{z}_{i}^{1},\mathbf{z}_{i}^{2})/\tau)+Neg}$
(3)

and $Neg$ is defined as

$Neg=\sum_{j\neq i}exp(sim(\mathbf{z}_{i}^{1},\mathbf{z}_{j}^{2})/\tau),$ (4)

where $sim(\cdot)$ is the similarity function and $\tau$ is temperature
scaling parameter.

##  III Proposed Method

Original HSI is 3-D structured and can be denoted as
$\mathcal{S}\in\mathbb{R}^{H\times{W}\times{C}}$ with resolution of
$H\times{W}$ and $C$ spectral bands. The target is to assign the total $N$
pixels into $K$ classes where $N=H\times{W}$. Our method consists of pixel-
level pre-training and superpixel-level contrastive clustering. Since
clustering is performed for superpixels, pixels’ labels are inherited from the
belonged superpixels. The overall architecture of the proposed method is shown
in Fig. 1. First, we employ a variational autoencoder composed of hybrid 3-D
and 2-D CNNs for pixel-level pre-training, extracting high-order spatial and
spectral features from HSI pixels. Subsequently, superpixel segmentation is
performed and the graph structure of superpixels is constructed. In superpixel
graph contrastive clustering, we propose two semantical-invariant data
augmentations for superpixels and employ sample-level alignment together with
clustering-center-level contrast as the target of network optimization.
Clustering results are obtained by K-means and then high-confidence clustering
centers are recomputed. We conduct clustering and network optimization
alternatively.

TABLE I: Architecture of Pixel Pre-training Network and -1 Denotes the Batch
Size.

Encoder |  |  | $Decoder$ |  |
---|---|---|---|---|---
Layer | Kernel | Output | $Layer$ | Kernel | Output
Input | - | [-1,1,30,27,27] | $Input$ | - | [-1,128]
Conv3D | [8,7,3,3] | [-1,8,24,25,25] | $FC$ | [128,256] | [-1,256]
BN3D | - | - | $FC$ | [256,23104] | [-1,23104]
Conv3D | [16,5,3,3] | [-1,16,20,23,23] | $Reshape$ | - | [-1,64,19,19]
BN3D | - | - | $DeConv2D$ | [576,3,3] | [-1,576,21,21]
Conv3D | [32,3,3,3] | [-1,32,18,21,21] | $BN2D$ | - | -
BN3D | - | - | $Reshape$ | - | [-1,32,18,21,21]
Reshape | - | [-1,576,21,21] | $DeConv3D$ | [16,3,3,3] | [-1,16,20,23,23]
Conv2D | [64,3,3] | [-1,64,19,19] | $BN3D$ | - | -
BN2D | - | - | $DeConv3D$ | [8,5,3,3] | [-1,8,24,25,25]
GAP | - | [-1,64,4,4] | $BN3D$ | - | -
Reshape | - | [-1,1024] | $DeConv3D$ | [1,7,3,3] | [-1,1,30,27,27]
FC | [1024,512] | [-1,512] | $BN3D$ | - | -
FC | [512,128] | [-1,128] |  |  |

###  III-A Pixel-level Pre-training

In this section, we perform pixel-level pre-training to fully extract spatial
and spectral features of HSI. Following the existing structure[35, 43], we
utilize 3-D and 2-D hybrid CNNs to construct variational autoencoder (VAE)[44]
and train the network without supervision. First, PCA is used to reduce the
number of bands for the reduction of computational cost. Then, we use a
sliding window to split the HSI into pixel cubes
$\mathbf{P}_{i}\in\mathbb{R}^{w\times{w}\times{h}}(i=1,2,...,N)$ where
$\mathbf{P}_{i}$ is the $i$-th pixel cube, $w$ and $h$ are the window size and
the number of reduced bands respectively. In this manner, pixels are
represented by a fixed-size neighborhood, which utilizes spatial information
to mitigate spectral variability. The encoder mainly consists of two layers of
3-D convolution and one layer of 2-D convolution. We add a global average
pooling layer (GAP) at the end to fix the output dimension. The output of GAP
for the $i$-th pixel cube is then reshaped into
$\mathbf{x}_{i}^{p}\in\mathbb{R}^{d}$. Next, we use fully-connected layers
(FC) to fit the mean vector $\bm{\mu}$ and standard deviation vector
$\bm{\sigma}$ of the data distribution, and the latent code $\mathbf{q}$ is
calculated with re-parameterization trick:

$\mathbf{q}=\bm{\mu}+\epsilon\times\bm{\sigma},$ (5)

where $\epsilon$ is sampled from standard normal distribution and
$\bm{\mu},\bm{\sigma},\mathbf{q}\in\mathbb{R}^{d^{\prime}}$. The decoder
consists of symmetric deconvolution layers to reconstruct the input pixel
cube, and its input and output are $\mathbf{q}$ and
$\overline{\mathbf{P}}_{i}$. In addition, we add batch normalization layers
after each convolution layer to reduce the risk of overfitting. The specific
structure of the pre-training network is shown in Table I.

Loss function of VAE is composed of distribution loss
$\mathcal{L}_{\bm{\mu},\bm{\sigma}^{2}}$ and reconstruction loss
$\mathcal{L}_{recon}$, which are calculated by (6) and (7):

$\mathcal{L}_{\bm{\mu},\bm{\sigma}^{2}}=\frac{1}{2}\sum_{i=1}^{d^{\prime}}(\mu_{(i)}^{2}+\sigma_{(i)}^{2}-log\sigma_{(i)}^{2}-1),$
(6)
$\mathcal{L}_{recon}=\frac{1}{2}\sum_{i=1}^{N}(\|\mathbf{P}_{i}-\overline{\mathbf{P}}_{i}\|_{F}^{2}),$
(7)

where the subscript $(i)$ denotes the $i$-th element in the vector, and the
total loss of pre-training is calculated by (8):

$\mathcal{L}_{pre-
train}=\mathcal{L}_{\bm{\mu},\bm{\sigma}^{2}}+\mathcal{L}_{recon}.$ (8)

After pre-training, we save the output of GAP
$\mathbf{X}^{p}=[\mathbf{x}_{1}^{p};\mathbf{x}_{2}^{p};...;\mathbf{x}_{N}^{p}]\in\mathbb{R}^{N\times
d}$ so that pixel’s representations with high-order spatial and spectral
information are accessible.

###  III-B Superpixel Segmentation and Graph Construction

The large number of pixels in HSI results in high computational complexity.
Benefiting from the superpixel segmentation technology, HSI can be divided
into several homogeneous regions. Specifically, we adopt entropy rate
superpixel segmentation (ERS) to divide HSI into $M$ non-overlapping
superpixels, denoted as $\mathcal{S}=\bigcup_{i=1}^{M}\mathbf{x}^{sp}_{i}$
where $\mathbf{x}^{sp}_{i}$ represents the $i$-th superpixel. As $M\ll N$, the
computational complexity can be greatly reduced by converting the clustering
task from the pixel level to the superpixel level.

To obtain superpixels’ representations, we first construct a map matrix
$\mathbf{T}\in\mathbb{R}^{M\times N}$ where $\mathbf{T}_{i,j}=1$ if the
$j$-pixel belongs to the $i$-th superpixel, otherwise 0.

Then each superpixel is represented as the average of pre-trained
representations of internal pixels by:

$\mathbf{X}^{sp}=\widetilde{\mathbf{T}}\mathbf{X}^{p},$ (9)

where $\mathbf{X}^{sp}\in\mathbb{R}^{M\times d}$ and $\widetilde{\mathbf{T}}$
is row-normalized $\mathbf{T}$. Similarity graph
$\mathbf{A}\in\mathbb{R}^{M\times M}$ is constructed based on the spatial
adjacent relationship between superpixels, i.e., $a_{m,n}=a_{n,m}=1$ if and
only if there exists at least one pair of
$(\mathbf{x}^{p}_{i},\mathbf{x}^{p}_{j})$ that is adjacent in HSI where
$\mathbf{x}^{p}_{i}$ belongs to the $m$-th superpixel and $\mathbf{x}^{p}_{j}$
belongs to the $n$-th superpixel.

For simplicity, we use $L$-layers graph convolutional network (GCN) as the
graph encoder to aggregate neighborhood information. The re-normalization
trick is applied by adding the self-loop:

$\widetilde{\mathbf{A}}=\mathbf{A}+\mathbf{I},$ (10)

where $\mathbf{I}$ is an identity matrix. The renormalized degree matrix
$\widetilde{\mathbf{D}}$ is calculated by
$\widetilde{d}_{ii}=\sum_{j}\widetilde{a}_{ij}$. Graph convolution is defined
as:

$\mathbf{H}_{(l)}=RELU(\widetilde{\mathbf{D}}^{-\frac{1}{2}}\widetilde{\mathbf{A}}\widetilde{\mathbf{D}}^{-\frac{1}{2}}\mathbf{H}_{(l-1)}\mathbf{W}_{(l)}),$
(11)

where $\mathbf{H}_{(l)}$ and $\mathbf{W}_{(l)}$ are the output and parameters
of the $l$-th layer and $\mathbf{H}_{(0)}=\mathbf{X}^{sp}$, $RELU(\cdot)$ is
non-linear activation function. We also denote the last layer’s output
$\mathbf{H}_{(L)}$ as $\mathbf{Z}$ for convenience.

###  III-C Superpixel Semantic-Invariant Data Augmentations

Here we propose two semantic-invariant augmentations for superpixels in
detail. First, we define pixel sampling augmentation (PSA). Each superpixel
contains several internal pixels, and its representation is the average of all
these pixels. It is a natural idea that pixels within a superpixel can serve
as semantic-invariant augmented views of this superpixel. Thanks to the high
precision of superpixel segmentation, internal pixels mostly belong to the
same class, and may not lead to semantic change (detailed information could be
found in Section IV-D). In practice, we randomly sample one pixel from each
superpixel as its positive augmentation view. All the sampled pixels are
denoted as $\mathbf{X}^{p}_{sampled}\in\mathbb{R}^{M\times d}$ which has the
same shape as $\mathbf{X}^{sp}$.

The second superpixel augmentation, model weight augmentation (MWA), is
enlightened by the multi-head idea[45, 46]. We design dual branches in the
last layer of the GCN encoder to obtain augmented views with semantic
consistency. Specifically, the $L$-th layer contains two branches with
unshared parameter matrices $\mathbf{W}_{(L)}^{1}$ and $\mathbf{W}_{(L)}^{2}$
which have the same architecture. Due to the different random initialization,
two branches’ outputs $\mathbf{H}_{(L)}^{1}$ and $\mathbf{H}_{(L)}^{2}$
represent two different views of the input superpixel but with the same
semantic, i.e., a further aggregation result based on the $(L-1)$-th layer.

Combining PSA with MWA, we can obtain different augmented views of the
superpixel. Specifically, representations of superpixels $\mathbf{X}^{sp}$ and
the corresponding sampled pixels $\mathbf{X}^{p}_{sampled}$ are the inputs of
the graph encoder. In the first $L-1$ layers, embeddings of superpixels and
sampled pixels are calculated by (11), denoted as $\mathbf{H}_{(L-1)}^{sp}$
and $\mathbf{H}_{(L-1)}^{p}$ respectively. In the last layer,
$\mathbf{H}_{(L-1)}^{sp}$ and $\mathbf{H}_{(L-1)}^{p}$ are inputted into both
two branches, thus the outputs include four embeddings:
$\mathbf{Z}^{sp_{1}},\mathbf{Z}^{sp_{2}},\mathbf{Z}^{p_{1}},\mathbf{Z}^{p_{2}}$
where $\mathbf{Z}^{sp_{1}}$ and $\mathbf{Z}^{p_{1}}$ are the outputs of the
first branch for superpixels and sampled pixels, while $\mathbf{Z}^{sp_{2}}$
and $\mathbf{Z}^{p_{2}}$ are of the second. Subsequently, we normalize each
row of these embeddings with $l_{2}$-norm.

Algorithm 1 SPGCC

0: HSI $\mathcal{S}$; Superpixel segmentation $\mathbf{T}$; Max iteration $I$.

0: Clustering results $\mathbf{R}$.

1: Split the input HSI $\mathcal{S}$ into pixel cubes $\mathbf{P}$ with a
sliding window.

2: Perform pixel-level pre-training on $\mathbf{P}$ to obtain pixel
representations $\mathbf{X}^{p}$.

3: Calculate the input superpixel features $\mathbf{X}^{sp}$ with (9) and get
adjacency matrix $\mathbf{A}$.

4: for $i=1$ to $I$ do

5: Randomly sample one pixel from each superpixel to get
$\mathbf{X}^{p}_{sampled}$.

6: Encode representations of superpixels and internal sampled pixels with (11)
to obtain
$\mathbf{Z}^{sp_{1}},\mathbf{Z}^{sp_{2}},\mathbf{Z}^{p_{1}},\mathbf{Z}^{p_{2}}$.

7: Concatenate $\mathbf{Z}^{sp_{1}}$ and $\mathbf{Z}^{sp_{2}}$ to get
$\mathbf{Z}^{sp}$.

8: Perform K-means on $\mathbf{Z}^{sp}$ to get clustering results $\mathbf{R}$
and centers $\mu_{k}$.

9: Select high-confidence samples according to (13).

10: Recompute high-confidence centers with (14).

11: Calculate sample-level alignment loss $\mathcal{L}_{SLA}$ with (18).

12: Calculate clustering-center-level contrast loss $\mathcal{L}_{CLC}$ with
(19).

13: Calculate the total loss $\mathcal{L}$ with (20).

14: Optimize model parameters with Adam.

15: end for

16: Map clustering results $\mathbf{R}$ to superpixel level to pixel level.

17: return $\mathbf{R}$

###  III-D Superpixel Contrastive Clustering

In this section, we jointly get clustering results and optimize superpixels’
embeddings under a contrastive learning framework. First, we concatenate the
two views $\mathbf{Z}^{sp_{1}}$ and $\mathbf{Z}^{sp_{2}}$ to get the final
embeddings of superpixels:

$\mathbf{Z}^{sp}=[\mathbf{Z}^{sp_{1}},\mathbf{Z}^{sp_{2}}].$ (12)

Then we perform K-means on $\mathbf{Z}^{sp}$ to obtain clustering assignments
and the $k$-th clustering center is denoted as $\bm{\mu}_{k}(k=1,2,...,K)$.
Considering the clustering centers are susceptible to outliers, to select
high-confidence samples and extract more reliable clustering information[46],
we calculate the distance between samples and their nearest clustering
centers:

$d_{i}=\min_{1\leq k\leq K}\|\mathbf{z}^{sp}_{i}-\bm{\mu}_{k}\|_{2}^{2},$ (13)

where $\mathbf{z}^{sp}_{i}$ and $d_{i}$ denote the $i$-th superpixel embedding
and its distance to the clustering center. Smaller distance means higher
clustering confidence, so superpixels are sorted in ascending order according
to their distances, and the top $\lambda$ are selected as high-confidence
samples. Indices of high-confidence samples in class $k$ are denoted as
$h_{k}=\\{h_{k}^{1},h_{k}^{2},...,h_{k}^{n_{k}}\\}$. Next, high-confidence
clustering centers of two views
$\mathbf{c}_{k}^{1},\mathbf{c}_{k}^{2}(k=1,2,...,K)$ that are more reliable
than $\bm{\mu}_{k}$ are recomputed only based on these high-confidence
samples:

$\mathbf{c}_{k}^{1}=\frac{1}{|h_{k}|}\sum_{i\in
h_{k}}\mathbf{z}_{i}^{sp_{1}},\mathbf{c}_{k}^{2}=\frac{1}{|h_{k}|}\sum_{i\in
h_{k}}\mathbf{z}_{i}^{sp_{2}}.$ (14)

Traditional contrastive clustering methods suffer from low intra-class
consistency caused by wrong negative pairs[47]. To overcome this deficiency,
we utilize sample-level alignment (SLA) and clustering-center-level contrast
(CLC) as optimization objectives. Briefly, at the sample level, we do not
treat different superpixels as negative pairs since they may belong to the
same class. Superpixel embeddings are only encouraged to be aligned with their
semantic-invariant augmented views. At the clustering center level, for a
certain clustering center, it is natural to set the other $K-1$ centers as its
negative views and perform contrast. Combining SLA with CLC, we hope to get
superpixel embeddings with better intra-class similarity and inter-class
dissimilarity.

For SLA, since we have obtained multiple augmented views of superpixels, we
perform six pairs of alignments in total, which can be divided into three
groups. The first group, $\mathbf{Z}^{sp_{1}}$ and $\mathbf{Z}^{sp_{2}}$,
$\mathbf{Z}^{p_{1}}$ and $\mathbf{Z}^{p_{2}}$, aims to encourage the outputs
of dual branches to be consistent. The second group, $\mathbf{Z}^{sp_{1}}$ and
$\mathbf{Z}^{p_{1}}$, $\mathbf{Z}^{sp_{2}}$ and $\mathbf{Z}^{p_{2}}$, aims to
pull the embeddings of superpixels close to sampled internal pixels. The third
group, $\mathbf{Z}^{sp_{1}}$ and $\mathbf{Z}^{p_{2}}$, $\mathbf{Z}^{sp_{2}}$
and $\mathbf{Z}^{p_{1}}$, combines the above two targets together. Loss
functions for three groups of alignments are defined as:

$\mathcal{L}_{SLA}^{group_{1}}=\|\mathbf{Z}^{sp_{1}}-\mathbf{Z}^{sp_{2}}\|_{F}^{2}+\|\mathbf{Z}^{p_{1}}-\mathbf{Z}^{p_{2}}\|_{F}^{2},$
(15)
$\mathcal{L}_{SLA}^{group_{2}}=\|\mathbf{Z}^{sp_{1}}-\mathbf{Z}^{p_{1}}\|_{F}^{2}+\|\mathbf{Z}^{sp_{2}}-\mathbf{Z}^{p_{2}}\|_{F}^{2},$
(16)
$\mathcal{L}_{SLA}^{group_{3}}=\|\mathbf{Z}^{sp_{1}}-\mathbf{Z}^{p_{2}}\|_{F}^{2}+\|\mathbf{Z}^{sp_{2}}-\mathbf{Z}^{p_{1}}\|_{F}^{2},$
(17)

and the total loss of SLA is:

$\mathcal{L}_{SLA}=\frac{1}{6}(\mathcal{L}_{SLA}^{group_{1}}+\mathcal{L}_{SLA}^{group_{2}}+\mathcal{L}_{SLA}^{group_{3}}).$
(18)

For CLC, we push away embeddings of different clustering centers, while
pulling together the same clustering center’s different views:

$\mathcal{L}_{CLC}=\frac{1}{K}\sum_{k=1}^{K}-log\frac{exp(\mathbf{c}_{k}^{1}\cdot\mathbf{c}_{k}^{2}/\tau)}{\sum_{j=1}^{K}exp(\mathbf{c}_{k}^{1}\cdot\mathbf{c}_{j}^{2}/\tau)},$
(19)

where $\tau$ is the temperature scaling parameter.

Finally, weight parameter $\alpha$ is adopted to balance the above two losses:

$\mathcal{L}=\mathcal{L}_{SLA}+\alpha\mathcal{L}_{CLC}.$ (20)

In this manner, clustering assignment and embedding optimization are performed
alternatively. The overall training detail is shown in Algorithm 1.

##  IV Results and Analysis

###  IV-A Datasets

We evaluate our method on three public HSI datasets, including India Pines,
Salinas, and Pavia University. An overview of three datasets can be found in
Table II.

TABLE II: Overview of India Pines, Salinas and Pavia University. Dataset |  Size | #Band | #Class | #Labeled pixels
---|---|---|---|---
India Pines |  (145, 145) | 200 | 16 | 10249
Salinas |  (512, 217) | 204 | 16 | 54129
Pavia University |  (610, 340) | 103 | 9 | 42776

1) India Pines: This dataset was captured by the AVIRIS sensor over North-
western Indiana in 1992. The size of this scene is 145 $\times$ 145 of which
10249 pixels are labeled. After removing 24 water absorption bands ([104-108],
[150-163], 220), each pixel is described by 200 remaining spectral bands. This
dataset contains 16 different land cover types.

2) Salinas: This dataset was captured by the AVIRIS sensor over Salinas
Valley, California. The resolution of this scene is 512 $\times$ 217 with
54129 labeled pixels. Similarly to India Pines, 20 water absorption bands
([108-112], [154-167], 224) are removed. Salinas contains 204 spectral bands
and 16 different land cover types.

3) Pavia University: This dataset was collected by the ROSIS sensor over
Pavia, northern Italy. It contains 103 spectral bands and 610 $\times$ 340
pixels, of which 42776 are labeled. In this scene, 9 different types of land
cover are recorded.

###  IV-B Experiment Setup

1) Implementation Details: Our method consists of pixel-level pre-training and
superpixel-level contrastive clustering. For the first stage, we divide 3-D
pixel cubes with a fixed window size of 27 $\times$ 27 and reduce the spectral
dimension to 30 on India Pines, 15 on Salinas, and 15 on Pavia University for
computational efficiency. The architecture of CNNs is the same as shown in
Table I for three datasets. Adam optimizer is utilized. We set the learning
rate and weight decay to $1\text{\times}{10}^{-3}$ and
$5\text{\times}{10}^{-4}$ respectively. When the pre-training of the VAE
module is finished, all the pixel cubes are mapped to feature vectors and
saved for the following learning stage. For the second stage, we use a
$3$-layers GCN as the graph encoder, where the sizes of the hidden layer and
the output layer are 1024 and 512 respectively. We perform K-means on
superpixel embeddings to obtain clustering assignments per five epochs. The
temperature scaling parameter of contrast is fixed to 0.5, and the weight of
clustering-center-level contrast loss is fixed to 0.1. We tune the learning
rate $\eta$ from $1\text{\times}{10}^{-5}$ to $1\text{\times}{10}^{-1}$. Other
hyperparameters are shown in Table III. All the experiments are repeated 10
times and we report the average value and standard deviation of clustering
metrics finally.

TABLE III: Hyperparameter Settings. Dataset |  $M$ | $\eta$ | $\alpha$ | $\lambda$
---|---|---|---|---
India Pines |  1100 | $1\text{\times}{10}^{-5}$ | 0.1 | 0.75
Salinas |  2700 | $1\text{\times}{10}^{-5}$ | 0.1 | 0.55
Pavia University |  2200 | $1\text{\times}{10}^{-4}$ | 0.1 | 0.25

2) Compared Methods: To comprehensively evaluate the performance and
applicability of the proposed SPGCC, we introduce ten baselines for
comparison, including K-means[12], fuzzy c-means (FCM)[13], sparse subspace
clustering (SSC)[14], scalable sparse subspace clustering (SSSC)[48],
selective sampling-based scalable sparse subspace clustering (S5C)[49],
superpixelwise PCA approach (SuperPCA)[50], graph convolutional optimal
transport (GCOT)[51], neighborhood contrastive subspace clustering (NCSC)[9],
superpixel-level global and local similarity graph-based clustering (SGLSC)[8]
and dual graph autoencoder (DGAE)[10]. Specifically, for K-means and FCM, we
perform clustering at the pixel level. For SSC, SSSC, S5C, clustering is
performed at the superpixel level to overcome the high memory cost of the
subspace clustering. For other methods, we follow the original settings
proposed in their papers. A brief summarization of the evaluated methods can
be found in Table IV.

TABLE IV: Properties of the proposed SPGCC and Compared Methods.

Method |  Superpixel | Subspace clustering | Graph learning | Deep learning
---|---|---|---|---
K-means |  |  |  |
FCM |  |  |  |
SSC |  $\checkmark$ | $\checkmark$ |  |
SuperPCA |  $\checkmark$ |  |  |
SSSC |  $\checkmark$ | $\checkmark$ |  |
S5C |  $\checkmark$ | $\checkmark$ |  |
GCOT |  | $\checkmark$ | $\checkmark$ |
SGLSC |  $\checkmark$ | $\checkmark$ | $\checkmark$ |
DGAE |  $\checkmark$ |  | $\checkmark$ | $\checkmark$
NCSC |  $\checkmark$ | $\checkmark$ |  | $\checkmark$
SPGCC |  $\checkmark$ |  | $\checkmark$ | $\checkmark$

3) Clustering Metrics: To quantitatively evaluate the clustering performance,
nine popular metrics are utilized: overall accuracy (OA), average accuracy
(AA), kappa coefficient (Kappa), normalized mutual information (NMI), adjusted
Rand index (ARI), F1-score (F1), Precision, Recall and Purity. Among these
metrics, OA, AA, NMI, F1, Precision, Recall, and Purity range in [0,1] while
Kappa and ARI range in [-1, 1]. For each metric, a higher score indicates a
better clustering result. Similar to the common process of clustering
evaluation, the Hungarian algorithm[52] is adopted to map the predicted labels
according to the ground truth before calculating metrics. In addition, we log
the time cost of each method.

TABLE V: Clustering Performance on HSI Datasets, India Pines, Salinas and
Pavia University are Abbreviated to IP, SA and PU Respectively, “OOM” Denotes
“Out of Memory”, and $\bullet$/$\circ$ Indicates whether SPGCC is
Statistically Superior/Inferior to the Compared Methods According to the
Pairwise $t$-Test at 0.05 Significance Level.

Dataset | Metric | K-means | $FCM$ | SSC | SuperPCA | SSSC | S5C | GCOT | SLGSC | DGAE | NCSC | SPGCC
---|---|---|---|---|---|---|---|---|---|---|---|---
IP | OA | .3737±.0050$\bullet$ |  $.3159\textpm.0035$∙ | .4769±.0060$\bullet$ | .5321±.0349$\bullet$ | .5130±.0238$\bullet$ | .4429±.0270$\bullet$ | .5030±.0010$\bullet$ | .5859±.0162$\bullet$ |  .5879±.0312$\bullet$ | .5342±.0185$\bullet$ | .6759±.0136
| AA | .4152±.0167$\bullet$ |  $.3541\textpm.0094$∙ | .4620±.0413$\bullet$ | .5099±.0724$\bullet$ | .5483±.0391$\bullet$ | .4876±.0300$\bullet$ | .3738±.0031$\bullet$ |  .5703±.0254$\bullet$ | .5218±.0471$\bullet$ | .4382±.0294$\bullet$ | .6274±.0346
| Kappa | .3045±.0045$\bullet$ |  $.2586\textpm.0035$∙ | .4173±.0079$\bullet$ | .4858±.0360$\bullet$ | .4551±.0264$\bullet$ | .3817±.0320$\bullet$ | .4245±.0011$\bullet$ |  .5486±.0170$\bullet$ | .5403±.0336$\bullet$ | .4769±.0195$\bullet$ | .6385±.0151
| NMI | .4394±.0034$\bullet$ |  $.4038\textpm.0013$∙ | .5738±.0071$\bullet$ |  .6639±.0139$\bullet$ | .5918±.0146$\bullet$ | .5558±.0114$\bullet$ | .5177±.0007$\bullet$ | .6395±.0065$\bullet$ | .6441±.0113$\bullet$ | .5675±.0104$\bullet$ | .6767±.0116
| ARI | .2186±.0018$\bullet$ |  $.1809\textpm.0013$∙ | .3136±.0063$\bullet$ | .4061±.0387$\bullet$ | .3270±.0163$\bullet$ | .3014±.0215$\bullet$ | .3322±.0004$\bullet$ | .4073±.0091$\bullet$ |  .4712±.0403$\bullet$ | .3910±.0101$\bullet$ | .5654±.0174
| F1 | .3068±.0039$\bullet$ |  $.2594\textpm.0012$∙ | .3923±.0064$\bullet$ | .4704±.0343$\bullet$ | .4052±.0126$\bullet$ | .3810±.0187$\bullet$ | .4291±.0003$\bullet$ | .4655±.0079$\bullet$ |  .5303±.0367$\bullet$ | .4624±.0093$\bullet$ | .6126±.0150
| Precision | .3340±.0055$\bullet$ |  $.3322\textpm.0013$∙ | .4200±.0082$\bullet$ | .5351±.0488$\bullet$ | .4288±.0265$\bullet$ | .4108±.0247$\bullet$ | .3634±.0004$\bullet$ |  .5819±.0168$\bullet$ | .5815±.0362$\bullet$ | .4842±.0099$\bullet$ | .6930±.0254
| Recall | .2840±.0106$\bullet$ |  $.2128\textpm.0011$∙ | .3683±.0118$\bullet$ | .4218±.0374$\bullet$ | .3850±.0122$\bullet$ | .3558±.0194$\bullet$ |  .5238±.0003$\bullet$ | .3881±.0077$\bullet$ | .4879±.0389$\bullet$ | .4426±.0128$\bullet$ | .5491±.0107
| Purity | .5117±.0021$\bullet$ |  $.4995\textpm.0033$∙ | .5705±.0119$\bullet$ | .6862±.0269$\bullet$ | .6031±.0336$\bullet$ | .5644±.0201$\bullet$ | .5471±.0009$\bullet$ | .6995±.0133$\bullet$ |  .7023±.0194$\bullet$ | .6109±.0175$\bullet$ | .7427±.0170
| Time (s) | 22.6 | $31.3$ | 20.3 | 6.4 | 2 | 1.4 | 146.1 | 108 | 25.5 | 77.5 | 31.1
SA | OA | .6697±.0013$\bullet$ |  $.5752\textpm.0312$∙ | .7188±.0092$\bullet$ | .7258±.0606$\bullet$ | .7587±.0170$\bullet$ | .6120±.0171$\bullet$ | OOM |  .8209±.0228$\bullet$ | .6683±.0267$\bullet$ | .7647±.0254$\bullet$ | .8447±.0102
| AA | .6856±.0055$\bullet$ |  $.6316\textpm.0270$∙ | .5794±.0113$\bullet$ | .6639±.0553$\bullet$ | .6330±.0252$\bullet$ | .6103±.0289$\bullet$ | OOM | .7125±.0071 | .6577±.0430$\bullet$ | .6109±.0093$\bullet$ | .7216±.0146
| Kappa | .6328±.0010$\bullet$ |  $.5363\textpm.0316$∙ | .6876±.0101$\bullet$ | .6936±.0692$\bullet$ | .7304±.0191$\bullet$ | .5765±.0186$\bullet$ | OOM |  .8011±.0241$\bullet$ | .6323±.0283$\bullet$ | .7374±.0272$\bullet$ | .8269±.0113
| NMI | .7245±.0015$\bullet$ |  $.6935\textpm.0099$∙ | .8051±.0065$\bullet$ | .8267±.0143$\bullet$ | .8344±.0118$\bullet$ | .7623±.0097$\bullet$ | OOM | .8757±.0058$\circ$ | .7803±.0151$\bullet$ | .8361±.0205$\bullet$ | .8688±.0056
| ARI | .5255±.0018$\bullet$ |  $.4507\textpm.0253$∙ | .6882±.0073$\bullet$ | .6639±.0647$\bullet$ | .6800±.0154$\bullet$ | .5014±.0135$\bullet$ | OOM |  .7164±.0227$\bullet$ | .5618±.0310$\bullet$ | .6902±.0561$\bullet$ | .7598±.0146
| F1 | .5762±.0020$\bullet$ |  $.5017\textpm.0249$∙ | .7214±.0064$\bullet$ | .7030±.0543$\bullet$ | .7149±.0135$\bullet$ | .5491±.0118$\bullet$ | OOM |  .7458±.0208$\bullet$ | .6069±.0314$\bullet$ | .7240±.0507$\bullet$ | .7850±.0130
| Precision | .5498±.0031$\bullet$ |  $.5500\textpm.0157$∙ | .6903±.0090$\bullet$ | .6311±.1043$\bullet$ | .6683±.0170$\bullet$ | .5852±.0183$\bullet$ | OOM |  .7330±.0239$\bullet$ | .5998±.0233$\bullet$ | .6719±.0344$\bullet$ | .7633±.0136
| Recall | .6053±.0086$\bullet$ |  $.4628\textpm.0401$∙ | .7555±.0071$\bullet$ | .8085±.0249 | .7689±.0155$\bullet$ | .5174±.0091$\bullet$ | OOM | .7604±.0374$\bullet$ | .6227±.0868$\bullet$ | .7860±.0742 | .8081±.0153
| Purity | .6736±.0018$\bullet$ |  $.6707\textpm.0096$∙ | .7613±.0142$\bullet$ | .7430±.0523$\bullet$ | .7715±.0152$\bullet$ | .7068±.0166$\bullet$ | OOM | .8357±.0134 | .7216±.0136$\bullet$ | .7779±.0240$\bullet$ | .8455±.0095
| Time (s) | 109.5 | $155.2$ | 11.7 | 92.1 | 12.4 | 2.4 | OOM | 47 | 360.2 | 2408.3 | 150.6
PU | OA | .5341±.0000$\bullet$ |  $.5162\textpm.0100$∙ | .5517±.0000$\bullet$ | .4666±.0177$\bullet$ | .4877±.0066$\bullet$ | .4601±.0058$\bullet$ | OOM | .5718±.0029$\bullet$ |  .5953±.0408$\bullet$ | OOM | .6268±.0048
| AA | .5259±.0001$\bullet$ |  $\uline{.5531\textpm.0085}$∙ | .4881±.0000$\bullet$ | .5205±.0228$\bullet$ | .5073±.0227$\bullet$ | .4606±.0084$\bullet$ | OOM | .5367±.0029$\bullet$ | .5180±.0303$\bullet$ | OOM | .5903±.0053
| Kappa | .4336±.0000$\bullet$ |  $.4209\textpm.0110$∙ | .4563±.0000$\bullet$ | .3610±.0187$\bullet$ | .4002±.0067$\bullet$ | .3476±.0036$\bullet$ | OOM | .4822±.0031$\bullet$ |  .4994±.0397$\bullet$ | OOM | .5491±.0056
| NMI | .5330±.0000$\bullet$ |  $.5284\textpm.0081$∙ | .5148±.0000$\bullet$ | .4424±.0004$\bullet$ | .4922±.0098$\bullet$ | .4741±.0028$\bullet$ | OOM |  .5389±.0043$\bullet$ | .5348±.0208$\bullet$ | OOM | .5792±.0114
| ARI | .3213±.0000$\bullet$ |  $.3052\textpm.0049$∙ | .3526±.0000$\bullet$ | .2542±.0100$\bullet$ | .2960±.0063$\bullet$ | .2897±.0127$\bullet$ | OOM | .4276±.0025$\bullet$ | .4483±.0515 | OOM | .4630±.0091
| F1 | .4543±.0000$\bullet$ |  $.4342\textpm.0065$∙ | .5004±.0000$\bullet$ | .3966±.0091$\bullet$ | .4239±.0045$\bullet$ | .4336±.0135$\bullet$ | OOM | .5408±.0016$\bullet$ | .5649±.0476 | OOM | .5679±.0070
| Precision | .5688±.0000$\bullet$ |  $.5737\textpm.0037$∙ | .5378±.0000$\bullet$ | .5092±.0079$\bullet$ | .5726±.0106$\bullet$ | .5249±.0050$\bullet$ | OOM |  .6711±.0042$\bullet$ | .6524±.0297$\bullet$ | OOM | .7132±.0116
| Recall | .3782±.0000$\bullet$ |  $.3495\textpm.0094$∙ | .4678±.0000$\bullet$ | .3249±.0090$\bullet$ | .3365±.0040$\bullet$ | .3695±.0166$\bullet$ | OOM | .4528±.0008$\bullet$ | .4995±.0557 | OOM | .4719±.0055
| Purity | .6961±.0000$\bullet$ |  $.6930\textpm.0095$∙ | .6387±.0000$\bullet$ | .5944±.0094$\bullet$ | .6593±.0095$\bullet$ | .6619±.0000$\bullet$ | OOM |  .7190±.0022$\bullet$ | .7057±.0163$\bullet$ | OOM | .7731±.0075
| Time (s) | 113.7 | $82.2$ | 16.4 | 87.4 | 7.7 | 2.9 | OOM | 195.6 | 202 | OOM | 103.5

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Figure 2: Clustering maps on India Pines (a) ground truth, (b) K-means, (c)
FCM, (d) SSC, (e) SuperPCA, (f) SSSC, (g) S5C, (h) GCOT, (i) SGLSC, (j) DGAE,
(k) NCSC, (l) SPGCC.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

Figure 3: Clustering maps on Salinas (a) ground truth, (b) K-means, (c) FCM,
(d) SSC, (e) SuperPCA, (f) SSSC, (g) S5C, (h) SGLSC, (i) DGAE, (j) NCSC, (k)
SPGCC.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

Figure 4: Clustering maps on Pavia University (a) ground truth, (b) K-means,
(c) FCM, (d) SSC, (e) SuperPCA, (f) SSSC, (g) S5C, (h) SGLSC, (i) DGAE, (j)
SPGCC.

###  IV-C Comparison with Other Methods

In Table V, we report the clustering performance and the running time of SPGCC
together with compared methods on India Pines, Salinas, and Pavia University.
For each clustering metric, the optimal value is marked in bold, and the
suboptimal value is underlined. Clustering map visualizations of different
methods are shown in Figs. 2, 3, 4. All the experiments were conducted in a
Linux server with an AMD EPYC 7352 CPU and an NVIDIA RTX 3090 GPU. CUDA was
available for deep learning methods except NCSC on Salinas because of the huge
cost of memory.

According to the quantitative metrics and clustering maps, it could be easily
observed that the proposed SPGCC method outperforms other compared methods on
three datasets significantly, especially on India Pines. For example, OA
obtained by SPGCC has an improvement of 8.8%, 2.4%, and 3.2% over the second-
best method on the three datasets, respectively. Moreover, all the nine
clustering metrics of SPGCC exceed those of the compared methods on India
Pines. Only NMI and Recall on Salinas and Recall on Pavia University do not
achieve the highest scores. Results of the significance test further validate
the superiority of SPGCC (i.e., all metrics on India Pines, five metrics on
Salinas, and six metrics on Pavia University have statistical superiority to
the suboptimal methods).

GCOT and NCSC fail to run on Salinas or Pavia University due to the extreme
cost of memory. GCOT does not use superpixel segmentation and therefore is not
suitable for large-scale HSI datasets, indicating the necessity of using
superpixel segmentation. Despite the risk of grouping pixels of different
classes into the same superpixel, superpixel segmentation brings high
homogeneity, and we have verified that the accuracy of segmentation is
generally above 99% (in Section IV-D). NCSC needs to perform convolution on
all pixel cubes in a single batch due to the network structure. Specifically,
although NCSC only conducts lightweight 2-D CNNs and introduces superpixel
segmentation to reduce the scale of subspace clustering, the pixel-level
convolutional network and the superpixel-level clustering network are not
separated, resulting in high computational costs. In contrast, the proposed
SPGCC also performs pixel-level convolution but in the pre-training stage
which is separated from the subsequent network. As a result, SPGCC could
extract high-order spectral and spatial features at the pixel level and
utilize superpixel segmentation to reduce data scale simultaneously.

Compared with the state-of-the-art graph learning methods (i.e., SGLSC and
DGAE), SPGCC performs better on all three datasets, with similar or less time
cost. SPGCC uses 3-D and 2-D hybrid CNNs at the pre-training stage to extract
high-order spatial and spectral features of HSI pixels, which is conducive to
subsequent clustering, while SGLSC performs clustering on raw HSI features.
Moreover, SPGCC is based on contrastive clustering and has an advantage over
the GAE-based method because of the clustering-oriented optimization target
(i.e., pulling superpixel representations of the same class close, while
pushing different clusters’ representations away). In this manner, SPGCC gains
clustering-friendly representations with better intra-class similarity and
inter-class dissimilarity.

(a)

(b)

(c)

Figure 5: The influence of different numbers of superpixels $M$ on (a) India
Pines, (b) Salinas, (c) Pavia University.

###  IV-D Parameters Analysis

Our method mainly contains four hyperparameters: number of superpixels ($M$),
learning rate ($\eta$), weight of clustering-center-level contrast loss
($\alpha$), and ratio of high-confidence samples ($\lambda$). They have been
tuned for different datasets. The adopted values are shown in Table III.

1) Impact of $M$: The number of superpixels $M$ controls the scale and
homogeneity of superpixels. The influence of $M$ on clustering performance is
shown in Fig. 5. It can be seen that the proposed SPGCC has robustness for
different $M$. As $M$ increases, different metrics almost show an increasing-
then-decreasing trend. The optimal $M$ are 1100, 2700, and 2200 for three
datasets, respectively.

In Fig. 6, we show the effect of different $M$ on the running time of our
algorithm. Settings of $M$ keep consistency with Fig. 5. Specifically, we fix
the training epochs to 200 and record the running time. Due to the efficient
parallelism of GCN, the increase in time consumption is lower than that of
$M$.

The impact of $M$ on superpixel segmentation accuracy $SP_{acc}$ is also
demonstrated in Fig. 6, which is defined as:

$SP_{acc}=\frac{\sum_{i=1}^{M}n_{i}^{dominant}}{\sum_{i=1}^{M}n_{i}^{labeled}},$
(21)

where $n_{i}^{dominant}$ and $n_{i}^{labeled}$ are the number of labeled
pixels in the dominant class and total labeled pixels within the $i$-th
superpixel, respectively. With the increase of $M$, $SP_{acc}$ shows a gradual
upward trend, indicating that pixels within one superpixel are more likely to
belong to the same class in a higher granularity partition. Since $SP_{acc}$
typically exceeds 99%, different $M$ have minimal impact on the accuracy of
superpixel segmentation.

(a)

(b)

(c)

Figure 6: Time cost and superpixel segmentation accuracy for different numbers
of superpixels $M$ on (a) India Pines, (b) Salinas, (c) Pavia University.

(a)

(b)

(c)

Figure 7: The influence of different ratio of selected high-confidence samples
$\lambda$ on (a) India Pines, (b) Salinas, (c) Pavia University.

(a)

(b)

(c)

Figure 8: The influence of different weight of clustering-center-level
contrast loss $\alpha$ on (a) India Pines, (b) Salinas, (c) Pavia University.

2) Impact of $\lambda$: $\lambda$ controls the ratio of selected samples when
recomputing high-confidence clustering centers. Fig. 7 shows the effect of
$\lambda$ on clustering metrics. According to Fig. 7, the proposed SPGCC is
quite robust to parameter $\lambda$. Specifically, fluctuations of the
clustering performance on the India Pines dataset are slightly greater than
those on Salinas and Pavia University. One possible reason is that the noise
intensity on India Pines is larger than that on other datasets and the number
of divided superpixels is smaller, so a larger $\lambda$ makes clustering
centers more susceptible to noise samples, while a smaller $\lambda$ means
only a small part of samples being selected which cannot fully represent
centers of certain classes. Tuning results show that the optimal values of
$\lambda$ for the three datasets are respectively 0.75, 0.55, and 0.25.

3) Impact of $\alpha$: $\alpha$ introduces the clustering-center-level
contrast loss. Fig. 8 suggests that as $\alpha$ increases from a very small
value, the improvement brought by the clustering-center-level contrast
gradually increases. When $\alpha$ is equal to 0.1, clustering performances on
three datasets are almost optimal, although it is slightly better when
$\alpha$ is set to 0.05 on Salinas. When $\alpha$ exceeds 0.1, clustering
performance declines significantly, indicating that a too-large weight for
clustering-center-level contrast loss is harmful to superpixel representation
learning. For the sake of simplicity, we finally fix $\alpha$ to 0.1.

TABLE VI: Ablation Study Results, OA is Repored.

Method |  India Pines | Salinas | Pavia University
---|---|---|---
w/o pixel-level pre-training |  .5354 | .6783 | .5049
w/o high-confidence centers |  .6548 | .8364 | .6218
w/o pixel sampling augmentation |  .6615 | .8419 | .6164
w/o model weight augmentation |  .6416 | .8247 | .5954
w/o sample-level alignment |  .6503 | .8418 | .5993
w/o clustering-center-level contrast |  .6534 | .8365 | .6137
SPGCC |  .6759 | .8447 | .6268

TABLE VII: Different Data Augmentations for Superpixels, OA is Repored, Pixel Sampling Augmentation and Model Weight Augmentation are Abbreviated to “PSA” and “MWA”, Bold and Underline Denote the Optimal and Suboptimal Values Respectively. Augmentation |  India Pines | Salinas | Pavia University
---|---|---|---
Node mask |  .6274 | .8225 | .5801
Node noise |  .6202 | .8217 | .5783
Node shuffling |  .6401 | .8230 | .5822
Edge permutation |  .6327 | .8377 | .5899
PSA |  .6416 | .8247 | .5954
MWA |  .6615 | .8419 | .6164
PSA+MWA |  .6759 | .8447 | .6268
Figure 9: t-SNE visualization on India Pines (first row), Salinas (second row)
and Pavia university (third row) by different features: raw features (first
column), pixel-level pre-trained features (second column) and SPGCC learned
features (third column).

###  IV-E Ablation Study

To validate the effectiveness of each module in the proposed SPGCC, including
pixel-level pre-training, recomputed high-confidence clustering centers, two
types of semantic-invariant superpixel augmentations (i.e., pixel sampling
augmentation and model weight augmentation), two clustering-oriented training
targets (i.e., sample-level alignment and clustering-center-level contrast),
we conduct a series of ablation experiments. According to Table VI, the
complete SPGCC achieves the best performance. In addition, we summarize the
following four observations. First, pixel-level pre-training is very important
and brings more than 10% performance improvement, indicating that high-order
spectral and spatial features extracted from pre-training are more conducive
to clustering. Second, recomputing clustering centers only with a certain
proportion of high-confidence samples can alleviate the influence of noise
samples and thus improve clustering performance. Third, using two types of
semantic-invariant superpixel data augmentations simultaneously performs
better than using only one, which is consistent with the current opinion in
self-supervised contrastive learning that multiple reasonable data
augmentations are beneficial to representation learning. Fourth, combining
sample-level alignment with clustering-center-level contrast performs better,
indicating that sample-level alignment helps to increase intra-class
similarity while clustering-center-level contrast helps to increase inter-
class dissimilarity, thus making representations easier to cluster.

###  IV-F Different Superpixel-level Data Augmentations

We also compared the proposed two semantic-invariant superpixel augmentations,
i.e., pixel sampling augmentation (PSA) and model weight augmentation (MWA),
with general graph data augmentation methods[41] which may hurt the node
semantic in Table VII. The compared data augmentations include node mask, node
noise, node shuffling, and edge permutation. We tune the changing rate from
10% to 30% and report their best performance.

According to the experimental results in Table VII, the best clustering result
is achieved by combining PSA and MWA. When only one type of augmentation is
used, MWA outperforms the others, while PSA ranks second on India Pines and
Pavia University, and third on Salinas, showing the advantage of proposed
semantic-invariant augmentations. In addition, general graph data augmentation
methods require additional hyperparameters to control the intensity of
augmentation, while PSA and MWA do not but still perform better, indicating
the convenience and effectiveness of the designed augmentations.

###  IV-G Visualization of Superpixel Representations with t-SNE

$t$-distributed stochastic neighbor embedding (t-SNE)[53] is used for
visualization of superpixel representations. We show the visualization results
of the raw HSI features, pixel-level pre-trained features, and embeddings
after SPGCC training in Fig. 9. The clustering performance of each type of
feature is reported as the best result among ten times K-means.

The raw features of superpixels are difficult to distinguish and are unevenly
distributed in feature space. The clustering performance of the pixel-level
pre-trained features has been improved, but since the target of the pre-
training stage is not clustering-oriented, the distribution of superpixel
representations is relatively loose. After SPGCC training, both intra-class
similarity and inter-class dissimilarity of the features have been improved.
It can be seen that the distribution of features in the same class is more
compact, while features from different classes have better discriminability.

##  V Conclusion

In this article, we proposed a superpixel graph contrastive clustering (SPGCC)
algorithm to learn clustering-friendly representations of HSI. In SPGCC, high-
order spatial and spectral information of HSI is extracted and preserved by
pixel-level pre-training with 3-D and 2-D hybrid CNNs. To conduct superpixel-
level contrastive clustering, we design two semantic-invariant superpixel data
augmentations, i.e., pixel sampling augmentation and model weight
augmentation, where internal pixels and the output of dual-branch GCN serve as
augmented views of superpixels. In addition, we recompute the high-confidence
center of each class to mitigate the influence of outliers. The optimization
target is composed of sample-level alignment and clustering-center-level
contrast to overcome the weakness of common contrastive learning that pushes
representations in the same class away. Experimental results demonstrated that
SPGCC outperforms compared methods on different HSI datasets significantly. In
the future, we hope to mine the more accurate correlation between superpixels
to further guide superpixel-level graph learning for HSI.

## References

  * [1] C. McCann, K. S. Repasky, R. Lawrence, and S. Powell, “Multi-temporal mesoscale hyperspectral data of mixed agricultural and grassland regions for anomaly detection,” _ISPRS J. Photogramm. Remote Sens._ , vol. 131, pp. 121–133, 2017.
  * [2] R. D. M. Scafutto, C. R. de Souza Filho, and W. J. de Oliveira, “Hyperspectral remote sensing detection of petroleum hydrocarbons in mixtures with mineral substrates: Implications for onshore exploration and monitoring,” _ISPRS J. Photogramm. Remote Sens._ , vol. 128, pp. 146–157, 2017.
  * [3] U. Heiden, W. Heldens, S. Roessner, K. Segl, T. Esch, and A. Mueller, “Urban structure type characterization using hyperspectral remote sensing and height information,” _Landscape Urban Plann._ , vol. 105, no. 4, pp. 361–375, 2012.
  * [4] J. Moraga and H. S. Duzgun, “Jigsawhsi: a network for hyperspectral image classification,” _arXiv:2206.02327_ , 2022.
  * [5] H. Liu, Y. Jia, J. Hou, and Q. Zhang, “Global-local balanced low-rank approximation of hyperspectral images for classification,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 32, no. 4, pp. 2013–2024, 2022.
  * [6] M. Li, Y. Liu, G. Xue, Y. Huang, and G. Yang, “Exploring the relationship between center and neighborhoods: Central vector oriented self-similarity network for hyperspectral image classification,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 33, no. 4, pp. 1979–1993, 2023.
  * [7] J. Fan, T. Chen, and S. Lu, “Superpixel guided deep-sparse-representation learning for hyperspectral image classification,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 28, no. 11, pp. 3163–3173, 2018.
  * [8] H. Zhao, F. Zhou, L. Bruzzone, R. Guan, and C. Yang, “Superpixel-level global and local similarity graph-based clustering for large hyperspectral images,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–16, 2022.
  * [9] Y. Cai, Z. Zhang, P. Ghamisi, Y. Ding, X. Liu, Z. Cai, and R. Gloaguen, “Superpixel contracted neighborhood contrastive subspace clustering network for hyperspectral images,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–13, 2022.
  * [10] Y. Zhang, Y. Wang, X. Chen, X. Jiang, and Y. Zhou, “Spectral–spatial feature extraction with dual graph autoencoder for hyperspectral image clustering,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 32, no. 12, pp. 8500–8511, 2022.
  * [11] Y. Ding, Z. Zhang, X. Zhao, Y. Cai, S. Li, B. Deng, and W. Cai, “Self-supervised locality preserving low-pass graph convolutional embedding for large-scale hyperspectral image clustering,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–16, 2022.
  * [12] T. Kanungo, D. Mount, N. Netanyahu, C. Piatko, R. Silverman, and A. Wu, “An efficient k-means clustering algorithm: analysis and implementation,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 24, no. 7, pp. 881–892, 2002.
  * [13] T.-n. Yang, C.-j. Lee, and S.-j. Yen, “Fuzzy objective functions for robust pattern recognition,” in _IEEE Int. Conf. Fuzzy Syst. (FUZZ-IEEE)_ , 2009, pp. 2057–2062.
  * [14] E. Elhamifar and R. Vidal, “Sparse subspace clustering: Algorithm, theory, and applications,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 35, no. 11, pp. 2765–2781, 2013.
  * [15] R. Vidal and P. Favaro, “Low rank subspace clustering (lrsc),” _Pattern Recognit. Lett._ , vol. 43, pp. 47–61, 2014.
  * [16] P. Ji, M. Salzmann, and H. Li, “Efficient dense subspace clustering,” in _IEEE Winter Conf. Appl. Comput. Vis. (WACV)_ , 2014, pp. 461–468.
  * [17] Y. Cai, Z. Zhang, Z. Cai, X. Liu, X. Jiang, and Q. Yan, “Graph convolutional subspace clustering: A robust subspace clustering framework for hyperspectral image,” _IEEE Trans. Geosci. Remote Sens._ , vol. 59, no. 5, pp. 4191–4202, 2021.
  * [18] N. Huang, L. Xiao, and Y. Xu, “Bipartite graph partition based coclustering with joint sparsity for hyperspectral images,” _IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens._ , vol. 12, no. 12, pp. 4698–4711, 2019.
  * [19] M.-Y. Liu, O. Tuzel, S. Ramalingam, and R. Chellappa, “Entropy rate superpixel segmentation,” in _Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2011, pp. 2097–2104.
  * [20] R. Achanta, A. Shaji, K. Smith, A. Lucchi, P. Fua, and S. Süsstrunk, “Slic superpixels compared to state-of-the-art superpixel methods,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 34, no. 11, pp. 2274–2282, 2012.
  * [21] H. Zhai, H. Zhang, P. Li, and L. Zhang, “Hyperspectral image clustering: Current achievements and future lines,” _IEEE Geosci. Remote Sens. Mag._ , vol. 9, no. 4, pp. 35–67, 2021.
  * [22] C. Hinojosa, E. Vera, and H. Arguello, “A fast and accurate similarity-constrained subspace clustering algorithm for hyperspectral image,” _IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens._ , vol. 14, pp. 10 773–10 783, 2021.
  * [23] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in _Proc. Int. Conf. Learn. Represent. (ICLR)_ , 2017.
  * [24] P. Velickovic, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph attention networks,” in _Proc. Int. Conf. Learn. Represent. (ICLR)_ , 2018.
  * [25] T. N. Kipf and M. Welling, “Variational graph auto-encoders,” in _Proc. Adv. Neural Inf. Process. Syst. (NeurIPS)_ , 2016.
  * [26] J. Xie, R. Girshick, and A. Farhadi, “Unsupervised deep embedding for clustering analysis,” in _Proc. 33rd Int. Conf. on Int. Conf. Mach. Learn. (ICML)_ , 2016, p. 478–487.
  * [27] Z. Peng, H. Liu, Y. Jia, and J. Hou, “Deep attention-guided graph clustering with dual self-supervision,” _IEEE Trans. Circuits Syst. Video Technol._ , vol. 33, no. 7, pp. 3296–3307, 2023.
  * [28] T. Chen, S. Kornblith, M. Norouzi, and G. Hinton, “A simple framework for contrastive learning of visual representations,” in _Proc. 33rd Int. Conf. on Int. Conf. Mach. Learn. (ICML)_ , 2020, p. 1597–1607.
  * [29] K. He, H. Fan, Y. Wu, S. Xie, and R. Girshick, “Momentum contrast for unsupervised visual representation learning,” in _Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2020, pp. 9726–9735.
  * [30] J. Grill, F. Strub, F. Altché, C. Tallec, P. H. Richemond, E. Buchatskaya, C. Doersch, B. Á. Pires, Z. Guo, M. G. Azar, B. Piot, K. Kavukcuoglu, R. Munos, and M. Valko, “Bootstrap your own latent - A new approach to self-supervised learning,” in _Proc. Adv. Neural Inf. Process. Syst. (NeurIPS)_ , 2020.
  * [31] Y. Li, P. Hu, Z. Liu, D. Peng, J. T. Zhou, and X. Peng, “Contrastive clustering,” in _Proc. AAAI Conf. Artif. Intell. (AAAI)_ , 2021, pp. 8547–8555.
  * [32] Y. Chen, H. Jiang, C. Li, X. Jia, and P. Ghamisi, “Deep feature extraction and classification of hyperspectral images based on convolutional neural networks,” _IEEE Trans. Geosci. Remote Sens._ , vol. 54, no. 10, pp. 6232–6251, 2016.
  * [33] Y. Li, H. Zhang, and Q. Shen, “Spectral–spatial classification of hyperspectral imagery with 3d convolutional neural network,” _Remote Sens._ , vol. 9, no. 1, 2017.
  * [34] M. Kanthi, T. H. Sarma, and C. S. Bindu, “A 3d-deep cnn based feature extraction and hyperspectral image classification,” in _IEEE India Geosci. Remote Sens. Symp. (InGARSS)_ , 2020, pp. 229–232.
  * [35] S. K. Roy, G. Krishna, S. R. Dubey, and B. B. Chaudhuri, “Hybridsn: Exploring 3-d–2-d cnn feature hierarchy for hyperspectral image classification,” _IEEE Geosci. Remote. Sens. Lett._ , vol. 17, no. 2, pp. 277–281, 2020.
  * [36] Y. Cai, Y. Liu, Z. Zhang, Z. Cai, and X. Liu, “Large-scale hyperspectral image clustering using contrastive learning,” in _Proc. ACM Int. Conf. Inf. Knowl. Manag. (CIKM)_ , vol. 3052, 2021.
  * [37] X. Chen and K. He, “Exploring simple siamese representation learning,” in _Proc. IEEE/CVF Conf. Comput. Vis. Pattern Recognit. (CVPR)_ , 2021, pp. 15 750–15 758.
  * [38] Y. You, T. Chen, Y. Sui, T. Chen, Z. Wang, and Y. Shen, “Graph contrastive learning with augmentations,” in _Proc. Adv. Neural Inf. Process. Syst. (NeurIPS)_ , 2020.
  * [39] Y. You, T. Chen, Y. Shen, and Z. Wang, “Graph contrastive learning automated,” in _Proc. Int. Conf. Mach. Learn. (ICML)_ , vol. 139, 2021, pp. 12 121–12 132.
  * [40] L. Wu, H. Lin, C. Tan, Z. Gao, and S. Z. Li, “Self-supervised learning on graphs: Contrastive, generative, or predictive,” _IEEE Trans. Knowl. Data Eng._ , vol. 35, no. 4, pp. 4216–4235, 2023.
  * [41] J. Zhou, C. Xie, Z. Wen, X. Zhao, and Q. Xuan, “Data augmentation on graphs: A technical survey,” _arXiv:2212.09970_ , 2022.
  * [42] P. Velickovic, W. Fedus, W. L. Hamilton, P. Liò, Y. Bengio, and R. D. Hjelm, “Deep graph infomax,” in _Proc. Int. Conf. Learn. Represent._ , 2019.
  * [43] Z. Cao, X. Li, Y. Feng, S. Chen, C. Xia, and L. Zhao, “Contrastnet: Unsupervised feature learning by autoencoder and prototypical contrastive learning for hyperspectral imagery classification,” _Neurocomputing_ , vol. 460, pp. 71–83, 2021.
  * [44] D. P. Kingma and M. Welling, “Auto-encoding variational bayes,” in _Proc. Int. Conf. Learn. Represent. (ICLR)_ , 2014.
  * [45] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, L. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Proc. Adv. Neural Inf. Process. Syst. (NeurIPS)_ , 2017, pp. 5998–6008.
  * [46] X. Yang, Y. Liu, S. Zhou, S. Wang, W. Tu, Q. Zheng, X. Liu, L. Fang, and E. Zhu, “Cluster-guided contrastive graph clustering network,” in _Proc. AAAI Conf. Artif. Intell. (AAAI)_ , 2023, pp. 10 834–10 842.
  * [47] Z. Huang, J. Chen, J. Zhang, and H. Shan, “Learning representation for clustering via prototype scattering and positive sampling,” _IEEE Trans. Pattern Anal. Mach. Intell._ , vol. 45, no. 6, pp. 7509–7524, 2023.
  * [48] X. Peng, H. Tang, L. Zhang, Z. Yi, and S. Xiao, “A unified framework for representation-based subspace clustering of out-of-sample and large-scale data,” _IEEE Trans. Neural Networks Learn. Syst._ , vol. 27, no. 12, pp. 2499–2512, 2016.
  * [49] S. Matsushima and M. Brbic, “Selective sampling-based scalable sparse subspace clustering,” in _Proc. Adv. Neural Inf. Process. Syst. (NeurIPS)_ , 2019, pp. 12 416–12 425.
  * [50] J. Jiang, J. Ma, C. Chen, Z. Wang, Z. Cai, and L. Wang, “Superpca: A superpixelwise pca approach for unsupervised feature extraction of hyperspectral imagery,” _IEEE Trans. Geosci. Remote Sens._ , vol. 56, no. 8, pp. 4581–4593, 2018.
  * [51] S. Liu and H. Wang, “Graph convolutional optimal transport for hyperspectral image spectral clustering,” _IEEE Trans. Geosci. Remote Sens._ , vol. 60, pp. 1–13, 2022.
  * [52] H. W. Kuhn, “The hungarian method for the assignment problem,” _Nav. Res. Logist._ , vol. 2, no. 1-2, pp. 83–97, 1955.
  * [53] L. van der Maaten and G. Hinton, “Visualizing data using t-sne,” _J. Mach. Learn. Res._ , vol. 9, no. 86, pp. 2579–2605, 2008.


# An impulsive geomagnetic effect from an early-impulsive flare

Hugh S. Hudson,1,2, Edward W. Cliver3, Lyndsay Fletcher1,4, Declan A. Diver1,
Peter T. Gallagher5, Ying Li6, Christopher M. J. Osborne1, Craig Stark1,Yang
Su6

1SUPA School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ
UK
2Space Sciences Laboratory, University of California, Berkeley, CA 94720 USA
3National Solar Observatory, 3665 Discovery Drive, Boulder CO 80303 USA
4Rosseland Centre for Solar Physics, University of Oslo, PO Box 1029 Blindern,
NO-0315 Oslo, Norway
5Astronomy & Astrophysics Section, DIAS Dunsink Observatory, Dublin Institute
for Advanced Studies, Dublin, D15 XR2R, Ireland
6Key Laboratory of Dark Matter and Space Astronomy, Purple Mountain
Observatory, Chinese Academy of Sciences,
Nanjing 210034, China<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

The geomagnetic “solar flare effect” (SFE) results from excess ionization in
the Earth’s ionosphere, famously first detected at the time of the Carrington
flare in 1859. This indirect detection of a flare constituted one of the first
cases of “multimessenger astronomy,” whereby solar ionizing radiation
stimulates ionospheric currents. Well-observed SFEs have few-minute time
scales and perturbations of $>10$ nT, with the greatest events reaching above
100 nT. In previously reported cases the SFE time profiles tend to resemble
those of solar soft X-ray emission, which ionizes the D-region; there is also
a less-well-studied contribution from Lyman $\alpha$. We report here a
specific case, from flare SOL2024-03-10 (M7.4), in which an impulsive SFE
deviated from this pattern. This flare contained an “early impulsive”
component of exceptionally hard radiation, extending up to $\gamma$-ray
energies above 1 MeV, distinctly before the bulk of the flare soft X-ray
emission. We can characterize the spectral distribution of this early-
impulsive component in detail, thanks to the modern extensive wavelength
coverage. A more typical gradual SFE occurred during the flare’s main phase.
We suggest that events of this type warrant exploration of the solar physics
in the “impulse response” limit of very short time scales.

###### keywords:

Sun: EUV – Sun: X-rays – Sun: corona

††pubyear: WWW††pagerange: An impulsive geomagnetic effect from an early-
impulsive flare–LABEL:lastpage

##  1 Introduction

The first reported optical observation of a flare (Carrington, 1859) coincided
with the first recorded “solar flare effect” (SFE) (Stewart, 1861). It
appeared in all three components of the geomagnetic field in close time
coincidence with the flare’s optical emission, andhad a magnitude of about 120
nT (of order 1% of the background geomagnetic field) mainly in the horizontal
component of the field. Since this time hundreds of SFEs have been detected,
as reviewed recently by Curto (2020). The SFE for the Carrington flare remains
topically important in view of the dangers posed by extreme space-weather
events (e.g., Cliver et al., 2022).

Any terrestrial magnetometer on the sunlit side of the Earth will detect an
SFE of sufficient magnitude to outweigh the complicated natural variability of
the geomagnetic field. The effects on the vector field at a given station
depend upon its location relative to the sub-solar point, as well as other
factors (Curto, 2020). At middle latitudes in the northern hemisphere, as with
the Carrington flare, the appearance of a typical SFE is a smooth depression
of the horizontal component generally following the time evolution of the
solar soft X-ray flux. Curto (2020) describes the complexity of the SFE
detection and morphology in greater detail.

The recent event SOL2024-03-10 (GOES class111This adheres to the current NOAA
calibration (e.g. Hudson et al., 2024). M7.4) displays the usual pattern of an
SFE, in that it has a slowly-developing ionospheric component roughly
coinciding with the soft X-rays. It is strikingly different, though, in that
it also shows an early sharp component in response to an “early impulsive”
hard X-ray burst. In early impulsive flares, the $>$25 keV hard X-ray flux
increase is delayed by less than 30 s after the soft X-ray onset (Sui et al.,
2007). Such flares can depart from the observed close relationship found in
the common Neupert effect (Neupert, 1968; Hudson, 1991; Dennis & Zarro, 1993).
The Neupert effect underpins standard flare models by linking the energy in
accelerated particles to that required for injecting hot plasma into coronal
flare loops.

The simplest description of flare time development consists of three phases:
onset, impulsive, and gradual, as recognized in soft X-ray emissions. The
onset phase (e.g., Hudson et al., 2021) typically is faint and not relevant to
SFE formation; the impulsive phase (Kane & Anderson, 1970) includes non-
thermal hard X-ray and microwave signatures, while the gradual phase consists
of hot plasma injected into the corona, and the well-known flare loops
(Neupert, 1968). The “early impulsive” character of the event described here
falls outside this standard pattern because it precedes the normal impulsive
development, which we identify with the Neupert effect. This distinctive
property is immediately evident in the hard X-ray spectrogram from the
ASO-S/HXI instrument (Zhang et al., 2019; Su et al., 2019) on board ASO-S (Gan
et al., 2023), as seen in Figure 1.

Figure 1: Background-subtracted ASO-S X-ray spectrogram of SOL2024-03-10,
showing the early-impulsive phase at 12:08-12:10 UT. The legend indicates the
HXI detectors used for the signal and for the background corrections. We
interpret the 10-20 keV X-ray response as the thermal soft X-rays of the
gradual phase, and the spikes at 30 keV over 12:10–12:13 UT as the normal
impulsive phase. The black line shows the GOES 1-8 Å flux.

In Section 2 we survey some of the databases containing information about the
early-impulsive epoch of SOL2024-03-10, describing a broad spectral energy
distribution (SED) covering some six decades from the UV to the MeV
$\gamma$-ray range. We note that at the high-energy end of the observed SED,
the photons penetrate the Earth’s atmosphere to heights well below the
ionospheric D-region. This layer absorbs most of the solar soft X-ray and
Lyman-$\alpha$ fluxes (e.g., Mitra, 1974). The $\gamma-ray$ photons cause
multiple ionizations via secondary processes. The “specific ionization”
approximation (e.g., Fermi, 1950) implies roughly 30 eV per ion pair; a single
MeV photon will thereby produce more than 104 secondary ions in the mesosphere
and down into the stratosphere at altitudes of 30-40 km. Moldavanov (2002)
describes the physics of a $\gamma$-ray component, which he terms a
“crotchet222Often spelled “crochet” (the French word for “hook”), an early
descriptive term for an SFE. impulse component” or CIC.

##  2 Observations

###  2.1 Overview of magnetic records

The subject flare SOL2024-03-10T12 (M7.4) occurred at S13W38 in NOAA AR 13599
and its SFE immediately marked it as unusual, as confirmed directly by the
hard X-ray spectrogram made available by the ASO-S database in near real time
(Figure 1). Many spacecraft and observatories detected various aspects of this
event. We use magnetometer data at 1-sec cadence from Dunsink Observatory,
nearly on the meridian of the subsolar point at about 20∘ to the S, and check
with 3-s data from Glasgow U̱niversity Observatory333www.magie.ie, with the
comparison shown in Figure 2. Our two geomagnetic datasets come from Dunsink
(N 53∘ 23′ 13.40′′, W 6∘ 20′ 15.17′′) and Glasgow (N 55∘ 51′ 26′′, W 4∘ 15′
7′′).

Figure 2: Comparison of magnetometer records from Dunsink (1 sec sampling) and
Glasgow (3 sec sampling) for SOL2024-03-10, showing residuals against daily
median values for the horizontal component $B_{hor}$. The dotted lines show
approximate maximum times for the fast SFE component (12:09:20 UT) and the
normal SFE (12:16 UT) respectively.

Generally this event had a stronger geomagnetic effect ($\Delta
B_{hor}\approx$ 7 nT) than other flares of this GOES class. A distinguishing
feature of this SFE is its fast component at about 12:09 UT, marked with the
left-hand dotted line in Figure 2. The data agree reasonably well and have
correlated irregular background variations at the nT level. The Glasgow 3-s
data are smoother than the Dunsink 1-s data, as would be expected from the
longer integrations. Similar two-step SFE time profiles were recorded at
Eskdalemuir and Hartland in the UK444https://imag-
data.bgs.ac.uk/GIN_V1/GINForms2..

###  2.2 Hodogram representation

During the SFE we can follow the disturbance vector via a hodogram, as shown
in Figure 3 here. This shows the angular motions of the total geomagnetic
vector. The two panels show the diurnal circulation of the Dunsink geomagnetic
vector and a blown-up region showing the time interval of the impulsive SFE
component. The diurnal motion of the direction of ${\bf B}$ shows a
counterclockwise pattern. The details of the SFE motions have interesting
differences: the impulsive component shows a tilt to the N at onset, followed
by a return on approximately the same path. Then, as shown by the colors in
the right panel of Figure 3, the field again exhibits a rotation, but roughly
in the opposite sense (clockwise) from that of the diurnal effect.

Figure 3: Hodogram plots of geomagnetic inclination (zenith = 90∘) against
azimuthal angle, measured from N to E for the Dunsink data. These are the
geomagnetic inclination and declination, respectively. The left panel shows
the full day, during which the field executes a complicated counter-clockwise
motion typical of the diurnal effect at this geographical point. The right
panel focuses on the flare interval so that one can identify the excess due to
the SFE. The colored time ranges separate 12:08 (red), 12:10 (orange), 12:16
(blue), up to 12:22 UT. Note the different deflection patterns for the early-
impulsive component (red) and the main SFE.

The hodogram suggests (from the Biot-Savart law) that the initial impulse in
${\bf B}$ results from a westward horizontal current if driven from the lower
latitudes nearer the subsolar point, which for this flare lies approximately
on the same meridian as the magnetic observatories.

###  2.3 Neupert Effect

The event importantly displays a clear distinction between the early-impulsive
component and the normal Neupert-effect component, as shown in the time-series
comparison in Figure 4 (left panel). This simplifies the relationship easily
seen in the ASO-S hard X-ray spectrogram (Figure 1), but now relating the
timing to the GOES/XRS soft channel (1-8 Å).

Figure 4: Left panel, time series of Fermi/GBM hard X-ray (31 keV) and GOES
soft X-ray (1-8 Å) fluxes up to the time of GOES maximum. The right panel
shows the correlation between the hard X-ray flux and the time derivative of
the GOES flux at 1-sec sampling, with the early-impulsive phase points in red
and the later points in blue. Both sets of points show strong correlations.

The right panel of Figure 4 directly compares the time derivative of soft
X-rays with the hard X-ray flux (Fermi/GBM NaI detector n0, energy channel
22). This derivative comparison is complementary to Neupert’s original
integral version, which related the time integral of impulsive microwave
emission to the soft X-ray time history. Both sets of points indeed correlate
well, with Pearson’s $\rho=[0.77,0.85]$ respectively for the early-impulsive
and impulsive phases. The clear distinction of the two sets of points implies
that the early-impulsive particle acceleration drives a significantly weaker
rate of mass transfer (evaporation) into the corona.

###  2.4 Spectral Energy Distribution

As seen, the isolated fast SFE corresponds to an unusual “early impulsive”
flare emission. In this section we describe the spectral energy distribution
via time series across a broad spectral range. Figure 5 shows representative
data, ordered top to bottom in decreasing photon energy:

  1. 1.

Fermi/GBM counting rates above 300 keV: the $\gamma$-ray continuum. This
extends up to at least 1.5 MeV.

  2. 2.

Fermi/GBM counting rates at 20 keV, the more typical hard X-ray continuum
defining the main impulsive phase.

  3. 3.

GOES XRS-B (1-8 Å, or roughly 2 keV): the hot thermal plasma in flare loops.

  4. 4.

SDO/EVE/ESP (18 nm, nominally 0.15 keV): the EUV line forest from flare loops
and footpoint regions.

  5. 5.

SDO/AIA 304 Å (nominally 40 eV): largely chromospheric.

  6. 6.

SDO/AIA 1600 and 1700 Å (nominally 7.5 eV): The corona/chromosphere interface
region (Simões et al., 2019).

  7. 7.

ASO-S/SDI Lyman-$\alpha$ (around 10 eV): largely chromospheric.

Note that coronal thermal emissions dominate the soft X-ray range (iii). The
line emissions in items (iv–vii) all show the impulsive character of the hard
X-rays, however, with time scales for flux variation $\tau=S/{\dot{S}}$ as
short as 5 s. The Lyman-$\alpha$ line is known to have both impulsive and
gradual signatures in other flare events (Milligan et al., 2020). Notably the
GOES/SXR data do not show rapid fluctuations directly because of the slowly-
cooling coronal hot plasma of the flare loop system (e.g., Acton et al.,
1982).

Figure 5: Data timeseries across a wide spectral range. The top two panels
show Fermi/GBM in counts; the GOES and ESP data are in irradiance units
(W/m2), while the lower three panels show excess fluxes normalized to the
event peak. The vertical dotted line shows a reference time for the early-
impulsive emissions.

The multiwavelength timeseries data for the early-impulsive phase (Figure 5)
reveal two important properties. First, all wavebands excepting the soft
X-rays from GOES have quite similar time profiles despite their different
emission mechanisms. Second, the Figure does not include a plot of visible
continuum (the white-light flare). We have examined the HMI pseudocontinuum,
which normally shows white-light flares well, and see only faint and diffuse
emissions with little energy and no hint of the fine spatial structures seen
in the EUV imaging (Section 2.6). The left panel of Figure 4 shows the X-ray
time development for the event, comparing normalized signals for one Fermi/GBM
NaI detector at about 31 keV, and the standard GOES-16/XRS 1-8 Å data. The
hard X-ray timing matches that of the fast SFE quite well, but the maximum
(negative) excursion of the slow SFE is delayed by about 3 minutes from the
1-8 Å suggesting that the main driver of this later component is at longer
wavelengths.

###  2.5 Time Variability

At the longest wavelengths the available data do not have sufficient time
resolution to follow the hard X-ray variability, but the ESP broadband EUV
time-series data (item iv in Figure 5) capture this with 1/4 s sampling
(Didkovsky et al., 2012). Accordingly, because of its relevance to
understanding causality in the solar sources of the radiation, we compared
Fermi 31-keV with the ESP 18-24 nm during the early impulsive phase in Figure
6. This gives ambiguous conclusions: of the five major peaks in hard X-rays,
the EUV appears to match to within about a second for the three peaks at about
12:08:45, 12:09:10, and 12:09:18 UT; but the ESP peak at 12:09:00 UT clearly
lags the hard X-rays by 2-3 s. The final hard X-ray peaks after about 12:09:25
have no detectable counterparts in the EUV at all.

Figure 6: Hard X-rays (31 keV at 1 s binning), blue, compared with EUV 18-24
nm at 1/4 s binning, red.

The early-impulsive SFE has the approximate duration of the multispectral
emissions, but may not have sufficient signal-to-noise ratio to show the fine
structure seen in the hard X-ray and EUV bands.

###  2.6 Flare imaging

The detection of the early-impulsive process at EUV wavelengths allows us to
locate the footpoint regions of the flaring structure. We sketch its relevant
properties with the overlay in Figure 7 here. This shows a cutout region of
the AIA 304 Å image at 12:09:20 UT (peak pixel brightness $\approx 2.6\times
10^{4}$ DN/sec), overlaid with an image at 12:14:05 UT (peak pixel brightness
$\approx 5\times 10^{3}$ DN/sec). The two panels show the two images, each
with the other as contours for ready comparison of the geometry. We list some
of the key features:

  * •

The early-impulsive image consists of compact brightenings that are basically
unresolved, and thus contain sub-arcsecond structure.

  * •

The early-impulsive brightenings occur in two regions, in a widespread double-
footpoint pattern.

  * •

These bright points are unresolved in time by the AIA image cadences.

  * •

The gradual-phase image contours lie interior to the early-impulsive
footpoints and develop a two-ribbon structure.

  * •

The prominent outer circular-ribbon structure evolves only slowly, not moving
significantly.

This complicated development contains many clues to the flaring structure. The
configuration suggests spine-and-fan connectivity for the magnetic field, with
the outer structure identifiable with the intersection of the dome and the
photosphere (Fletcher et al., 2001; Masson et al., 2009; Reid et al., 2012).
In this interpretation, the images suggest that the main flare development
lies inside (below) the preflare dome, while the western early-impulsive
footpoint could show the location of the interior spine.

Figure 7: AIA 304 Å images ($64^{\prime\prime}\times 64^{\prime\prime}$) of
SOL2024-03-10 at 12:09:20 UT, just at the end of the early-impulsive peak
(color), and a main-phase image at 12:14:05 UT. The images are reversed-color
and independently scaled by the square root of the brightness. Left, the
early-impulsive image, showing compact footpoints, with contours from the
main-phase image; right, the reverse view. The large-scale arc, which we
suggest corresponds to the dome of the spine-and-fan geometry, does not move
greatly as the flare develops, and the main-phase emission concentrates
between the early-impulsive footpoints (see text for commentary).

##  3 Discussion

The early-impulsive flare SOL2024-03-10 gives us a chance to characterize an
unusual perturbation of the Earth’s upper atmosphere via its observed spectral
energy distribution. This suggests contributions across a broad altitude range
in the upper stratosphere and mesosphere, and this “test pulse” distinctly
differs from the usual SFE behavior both in the time domain and in the
spectrum. It should therefore make an interesting exercise for model
development of the current system (Moldavanov, 2002). We note that the
hodogram analysis reveals a different pattern than for the subsequent main
phase of the SFE, possibly implying some feedback in the development of the
enhanced currents.

Similarly, isolating the early-impulsive solar emission, with its relatively
hard spectrum, should interest modelers of flare development The relative
weakness of evaporation for a harder primary electron spectrum results from
the greater depth of penetration of higher-energy electrons, which can be
modeled by 1D radiation hydrodynamics (e.g. Reep et al., 2015). The major
puzzle presented by this event is its lack of white-light continuum, which
conventionally should result from excitation in the deep solar atmosphere
(e.g., Martínez Oliveros et al., 2012).

The impulsive SFE appears most strongly in the geomagnetic X component (N
pointing). According to the Biot-Savart law, this corresponds to a westward
horizontal current, consistent with a flare subsolar point on the meridian
(e.g. Curto, 2020). In the subsequent development of the normal SFE component
here, the hodogram shows a similar counter-clockwise circulation as seen in
the diurnal motion at Dunsink and Glasgow on the day of this event. It is
beyond the scope of this paper to study this SFE in detail, but many
geomagnetic observatories will have records of it at high time resolution.

##  4 Conclusions

We have reported on the geomagnetic effects of a recent flare exhibiting a
rapid ionospheric response to an “early impulsive” flare emission. The
spectral energy distribution for the solar emissions extend to the
$\gamma$-ray energy range. There are other such two-step SFE events,
presumably with to early $\gamma$-ray emission, such as SOL1984-04-24 (Cliver
et al., 2020), and so the mechanism isolated in our case is likely not unique.
Future modeling of SFE development should benefit from this new kind of test
pulse, with its reduced confusion from the usual soft X-ray component.

“Early impulsive” behavior is not rare (Sui et al., 2007), but we have shown
that such a burst may result in a deviation from the simplest form of the
Neupert effect (e.g. Veronig et al., 2005). In this case we find that the
early-impulsive emissions appear in compact footpoint regions, presumably
connected by large-scale fields in a spine-and-fan geometry.

##  5 Acknowledgements:

We thank the anonymous reviewer for helpful comments. The recognition of the
event we have studied was facilitated by the availability of ASO-S near-real-
time data via http://sprg.ssl.berkeley.edu/~tohban/browser/, for which we
thank Albert Shih. ASO-S is supported by the National Key R&D Program of China
2022YFF0503002, the Strategic Priority Research Program of the Chinese Academy
of Sciences, Grant No. XDB0560000, the NSFC 12333010 and NSFC 11820101002. L.
F. acknowledges support from grant ST/X000990/1 made by the UK’s Science and
Technical Facilities Council. H. S. H. thanks the University of Glasgow for
hospitalilty, and Graham Woan for great assistance with the Glasgow
magnetometer. For the purpose of open access, the authors have applied a
Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript
version arising from this submission.

##  6 Data Availability Statement:

All data used in this article are publicly available. We thank the
Magnetometer Network of Ireland (https://www.magie.ie) for the geomagnetic
record used here.

## References

  * Acton et al. (1982) Acton L. W., Leibacher J. W., Canfield R. C., Gunkler T. A., Hudson H. S., Kiplinger A. L., 1982, ApJ, 263, 409
  * Carrington (1859) Carrington R. C., 1859, MNRAS, 20, 13
  * Cliver et al. (2020) Cliver E. W., Hayakawa H., Love J. J., Neidig D. F., 2020, ApJ, 903, 41
  * Cliver et al. (2022) Cliver E. W., Schrijver C. J., Shibata K., Usoskin I. G., 2022, Living Reviews in Solar Physics, 19, 2
  * Curto (2020) Curto J. J., 2020, Journal of Space Weather and Space Climate, 10, 27
  * Dennis & Zarro (1993) Dennis B. R., Zarro D. M., 1993, Sol. Phys., 146, 177
  * Didkovsky et al. (2012) Didkovsky L., Judge D., Wieman S., Woods T., Jones A., 2012, Sol. Phys., 275, 179
  * Fermi (1950) Fermi E., ed. 1950, Nuclear physics
  * Fletcher et al. (2001) Fletcher L., Metcalf T. R., Alexander D., Brown D. S., Ryder L. A., 2001, ApJ, 554, 451
  * Gan et al. (2023) Gan W., et al., 2023, Sol. Phys., 298, 68
  * Hudson (1991) Hudson H. S., 1991, in Bulletin of the American Astronomical Society. p. 1064
  * Hudson et al. (2021) Hudson H. S., Simões P. J. A., Fletcher L., Hayes L. A., Hannah I. G., 2021, MNRAS, 501, 1273
  * Hudson et al. (2024) Hudson H., Cliver E., White S., Machol J., Peck C., Tolbert K., Viereck R., Zarro D., 2024, Sol. Phys., 299, 39
  * Kane & Anderson (1970) Kane S. R., Anderson K. A., 1970, ApJ, 162, 1003
  * Martínez Oliveros et al. (2012) Martínez Oliveros J.-C., et al., 2012, ApJ, 753, L26
  * Masson et al. (2009) Masson S., Pariat E., Aulanier G., Schrijver C. J., 2009, ApJ, 700, 559
  * Milligan et al. (2020) Milligan R. O., Hudson H. S., Chamberlin P. C., Hannah I. G., Hayes L. A., 2020, Space Weather, 18, e02331
  * Mitra (1974) Mitra A. P., 1974, Ionospheric effects of solar flares. Vol. 46, doi:10.1007/978-94-010-2231-6,
  * Moldavanov (2002) Moldavanov A. V., 2002, Journal of Physics D Applied Physics, 35, 1311
  * Neupert (1968) Neupert W. M., 1968, ApJ, 153, L59
  * Reep et al. (2015) Reep J. W., Bradshaw S. J., Alexander D., 2015, ApJ, 808, 177
  * Reid et al. (2012) Reid H. A. S., Vilmer N., Aulanier G., Pariat E., 2012, A&A, 547, A52
  * Simões et al. (2019) Simões P. J. A., Reid H. A. S., Milligan R. O., Fletcher L., 2019, ApJ, 870, 114
  * Stewart (1861) Stewart B., 1861, Philosophical Transactions of the Royal Society of London Series I, 151, 423
  * Su et al. (2019) Su Y., et al., 2019, Research in Astronomy and Astrophysics, 19, 163
  * Sui et al. (2007) Sui L., Holman G. D., Dennis B. R., 2007, ApJ, 670, 862
  * Veronig et al. (2005) Veronig A. M., Brown J. C., Dennis B. R., Schwartz R. A., Sui L., Tolbert A. K., 2005, ApJ, 621, 482
  * Zhang et al. (2019) Zhang Z., et al., 2019, Research in Astronomy and Astrophysics, 19, 160