Unsupervised Segmentation of Hyperspectral Images Using 3D Convolutional Autoencoders

TL;DR

This paper presents an unsupervised deep learning method using 3D convolutional autoencoders combined with clustering to segment hyperspectral images effectively without requiring labeled data.

Contribution

The authors introduce a novel unsupervised deep architecture that couples 3D convolutional autoencoders with clustering for hyperspectral image segmentation.

Findings

Achieves high-quality segmentation without prior class labels

Performs well on benchmark and real-life datasets

Outperforms traditional methods in unsupervised segmentation

Abstract

Hyperspectral image analysis has become an important topic widely researched by the remote sensing community. Classification and segmentation of such imagery help understand the underlying materials within a scanned scene, since hyperspectral images convey a detailed information captured in a number of spectral bands. Although deep learning has established the state of the art in the field, it still remains challenging to train well-generalizing models due to the lack of ground-truth data. In this letter, we tackle this problem and propose an end-to-end approach to segment hyperspectral images in a fully unsupervised way. We introduce a new deep architecture which couples 3D convolutional autoencoders with clustering. Our multi-faceted experimental study---performed over benchmark and real-life data---revealed that our approach delivers high-quality segmentation without any prior class labels.

Figures

Tables

Table 1. TABLE I: The results obtained over all benchmarks.

Unsupervised segmentation measures
Set$\to$	Sa		IP		PU
Algorithm$\downarrow$	NMI	ARS	NMI	ARS	NMI	ARS
1D-CNN*	0.885	0.725	0.705	0.586	0.786	0.771
1D-CNN* with PCA	0.860	0.685	0.640	0.493	0.360	0.215
1D-CNN* with ICA	0.873	0.723	0.641	0.502	0.616	0.556
1D-CNN* (S-MSI)	0.880	0.723	0.718	0.609	0.784	0.759
1D-CNN* (CAE)	0.886	0.738	0.663	0.496	0.748	0.658
GM	0.819	0.642	0.445	0.229	0.514	0.290
GM with PCA	0.830	0.654	0.443	0.235	0.530	0.404
GM with ICA	0.838	0.665	0.436	0.212	0.522	0.396
GM (S-MSI)	0.848	0.673	0.456	0.248	0.532	0.407
GM (CAE)	0.628	0.475	0.435	0.289	0.480	0.459
\hdashline$k$-means	0.732	0.538	0.437	0.211	0.546	0.350
$k$-means with PCA	0.724	0.524	0.430	0.204	0.545	0.324
$k$-means with ICA	0.730	0.535	0.381	0.178	0.477	0.263
$k$-means (S-MSI)	0.712	0.496	0.430	0.208	0.546	0.325
$k$-means (CAE)	0.710	0.503	0.451	0.297	0.539	0.336
\hdashline3D-CAE	0.714	0.533	0.431	0.231	0.553	0.339
3D-CAE with PCA	0.746	0.527	0.467	0.263	0.639	0.546
3D-CAE with ICA	0.839	0.644	0.504	0.278	0.538	0.316
3D-CAE (S-MSI)	0.728	0.531	0.442	0.241	0.601	0.450

Table 2. TABLE II: The ranking (averaged across all benchmarks), and the execution time of all methods. The best ranking is boldfaced.

	Ranking		Time (min)
Algorithm$\downarrow$	NMI	ARS	Sa	IP	PU	Mu
GM	7.33	9.00	11.83	1.63	2.89	91.49
GM with PCA	6.67	5.00	1.35	0.09	1.65	32.97
GM with ICA	7.67	6.00	0.87	0.29	0.50	5.40
GM (S-MSI)	4.33	3.33	1.03	0.15	1.41	21.33
GM (CAE)	12.33	6.00	0.98	0.18	0.92	18.42
\hdashline$k$-means	6.50	8.00	1.12	0.18	0.79	52.01
$k$-means with PCA	9.50	11.67	0.28	0.07	0.37	16.96
$k$-means with ICA	12.00	11.67	0.35	0.04	0.99	13.84
$k$-means (S-MSI)	9.67	11.67	0.31	0.06	0.34	19.75
$k$-means (CAE)	8.00	7.33	0.26	0.08	0.44	19.82
\hdashline3D-CAE	8.33	8.00	91.31	14.75	102.12	1994.20
3D-CAE with PCA	3.00	5.00	20.15	7.18	36.28	640.45
3D-CAE with ICA	3.67	6.33	16.88	5.37	27.99	573.52
3D-CAE (S-MSI)	6.00	6.00	20.33	5.36	50.60	553.55

Equations

$$L_r=\frac{1}{p}\sum _{i=1}^p{\left|\left|{\mathbf{x}}_i-{{\mathbf{x}}^{\prime }}_i\right|\right|}_2^2,$$

$$q_{ij}=\frac{(1+{\left|\left|z_i-{\mu }_j\right|\right|}^2{)}^{-1}}{{\sum }_{j=1}^J\left(1+{\left|\left|z_i-{\mu }_j\right|\right|}^2\right)}.$$

$$L_c=\mathrm{KL}(\mathcal{T}||q)=\sum _{i=1}^{p^{\prime }}\sum _{j=1}^J{t}_{ij}\log \frac{t_{ij}}{q_{ij}},$$

$$t_{ij}=\frac{q_{ij}^2/{\sum }_{i=1}^{p^{\prime }}{q}_{ij}}{{\sum }_{j=1}^J\left(q_{ij}^2/{\sum }_{i=1}^{p^{\prime }}{q}_{ij}\right)}.$$

$$L=L_r+\alpha L_c,$$

$$\mathrm{NMI}=\frac{\mathrm{MI}(A,B)}{\left[H(A)+H(B)\right]/2},$$

$$\mathrm{ARS}=\frac{\left(\genfrac{}{}{0ex}{}{n}{2}\right)\left(a+d\right)-\left[\left(a+b\right)\cdot \left(a+c\right)+\left(c+d\right)\cdot \left(b+d\right)\right]}{{\left(\genfrac{}{}{0ex}{}{n}{2}\right)}^2-\left[\left(a+b\right)\cdot \left(a+c\right)+\left(c+d\right)\cdot \left(b+d\right)\right]},$$

Full text

I Introduction

Hyperspectral imaging (HSI) provides detailed information about the material within a captured scene. It registers a number of spectral bands, commonly up to hundreds of them, and can be exploited to understand the location and characteristics of the objects in the process of HSI classification and segmentation. In classification, we assign a single label to an input HSI pixel, whereas in segmentation we are focused on finding the boundaries of objects within an image111Therefore, segmentation involves classification of separate pixels.. Due to the increased availability of hyperspectral sensors, HSI analysis has become an important research topic tackled by the machine learning, remote sensing, and pattern recognition communities. Such imagery has multiple applications in a plethora of fields, including biochemisty, biology, medicine, geosciences, military defense, food quality management and monitoring, pharmacy, and many more [1]. Hyperspectral imaging is an indispensable tool in Earth observation, as it captures Earth peculiarities that are useful in precision agriculture, managing environmental disasters, military defense applications, soil monitoring or prediction of environmental events.

HSI classification and segmentation techniques can be divided into conventional machine learning algorithms, requiring feature engineering (i.e., feature extraction and selection) [2, 3], and deep learning-powered approaches, in which the appropriate representation is learned during the training [4, 5, 6, 7, 8, 9, 10]. The success of deep learning is reflected in a variety of fields, where it established the state of the art. However, to deploy deep models in practice, we need large and representative ground-truth training sets. It is a serious limiting factor in hyperspectral Earth observation analysis, where transferring HSI from an imaging satellite back to Earth is extremely costly. Creating new ground-truth datasets is error-prone and requires building a thorough understanding of the materials within a scene. Therefore, it involves acquiring observational ground-sensor data—it is often cost- and time-inefficient. These difficulties result in a very small number of ground-truth HSIs. In [11], we analyzed 17 recent papers in which only seven benchmarks were exploited, with only three of them being “widely-used”: Pavia University (15 papers), Indian Pines (in 8 papers), and Salinas Valley (5 papers).

There are three main approaches to deal with the limited ground-truth hyperspectral sets: (i) data augmentation, (ii) transfer learning, and (iii) unsupervised analysis of HSI, with (i) and (ii) being exploited mostly in deep learning-powered techniques. Data augmentation is a process of generating artificial examples following the original data distribution. Such samples can extend the training sets or they can be elaborated at the inference time, to build an intrinsic ensemble-like deep model [12]. In transfer learning, we train feature extractors over the source training data, and apply it to the target data [13]. This approach allows us to benefit from the available data to train efficient extractors—the classification part of a network is later fine-tuned over the target (much smaller) HSI. Both augmentation and transfer learning require the annotated target sets to either input them to an augmentation engine, or to utilize them for fine-tuning the deep models. Hence, their usefulness is limited in scenarios where manually-analyzed HSIs do not exist and are infeasible to generate.

On the other hand, unsupervised segmentation offers the possibility of processing HSI without any prior class labels. Although the literature in unsupervised HSI segmentation is rather limited, there exist conventional machine learning approaches (exploiting hand-crafted features) which benefit from mean shift filtering [14], diffusion-based dimensionality reduction followed by clustering [15], and phase-correlation analysis[16]. In [17], the authors used a fully conv-deconv deep network for unsupervised spectral-spatial feature learning. Then, the convolutional subnetwork was used as a generic feature extractor over the target (labeled) data in the transfer learning approach. A similar technique of extracting deep features using stacked sparse autoencoders, and later embedding them into linear support vector machines has been proposed in [18].

In this letter, we tackle the problem of limited ground-truth hyperspectral sets, and propose a deep learning technique for unsupervised HSI segmentation. Inspired by a recent work by Guo et al. [19], we introduce a 3D convolutional autoencoder architecture to learn embedded features, which later undergo clustering (Section II). This clustering is performed during the network training with a clustering-oriented loss, therefore our method delivers end-to-end unsupervised HSI segmentation. To the best of our knowledge, such approaches have not been investigated in the HSI literature so far. We performed a multi-faceted experimental study—over benchmark and real-life hyperspectral data—to understand the abilities of the proposed technique. It showed that our method offers high-quality and consistent segmentation, and does not require any prior class labels to effectively segment HSI (Section III).

II Unsupervised HSI Segmentation Using 3D Convolutional Autoencoders

In our unsupervised HSI segmentation approach (3D-CAE) inspired222In [19], the authors used much deeper architectures without dropout over full input images (with 2D kernels only) in the context of image classification. by [19], we exploit 3D convolutional autoencoders (CAEe) to extract deep features which later undergo clustering (Fig. 1). In the encoding part of the network, we capture both spectral and spatial features within an input three-dimensional hyperspectral patch $\bm{x}$ (of size $5\times 5\times{\rm B}$, where ${\rm B}$ is the number of bands; the patches are extracted with unit stride) using two convolutional layers denoted as ${\rm Conv}@h_{k}\times w_{k}\times d_{k}\times k$, for which we define the height ($h_{k}$), width ($w_{k}$), and depth ($d_{k}$) of the kernels, alongside the number of kernels in this layer (here, $k=32$ for all convolution/transposed convolution layers)—each kernel moves with unit stride in each direction. These convolutional layers are interleaved with one dropout layer (with the dropout probability of $p=0.5$) acting as a regularizer. The central-pixel features in the patch are later re-shaped to form a 1D vector which becomes an input to a fully-connected (embedding) layer with $n=25$ neurons, whose output is the latent vector. These embedded features are transformed back to the original 3D patch (to get the output 3D patch $\bm{x^{\prime}}$) in the decoding part of a CAE which is a mirrored version of the encoder with the transposed convolutions applied for up-sampling. The CAE is learned in the first training stage with the following reconstruction loss:

[TABLE]

where $p$ is the number of 3D patches in a batch. This stage, in which we do not use the clustering layer, runs until reaching convergence or the stopping condition (in this letter, the optimization terminates if the difference between two consecutive reconstruction loss values is less than $\epsilon=10^{-6}$).

In the second training stage, we modify the loss function and “switch on” the clustering layer—it is connected to the embedded layer of CAE which outputs the latent vector $z_{i}$ for the $i$-th input patch. The embedded features are being assigned a soft label $q_{i}$ in the clustering layer. As proposed in [19], this layer maintains the cluster centers $\mu_{j}$, where $j=1,2,\dots J$, and $J$ is the number of clusters, as trainable weights. The probability of assigning an input 3D patch $\bm{x}_{i}$ to each $j$-th cluster is generated using the Student’s $t$-distribution:

[TABLE]

Finally, the clustering loss is:

[TABLE]

where $p^{\prime}$ is the number of pixels in the batches, KL is the Kullback-Leibler divergence, and $\mathcal{T}$ is the target distribution:

[TABLE]

The clustering loss is incorporated into the total loss function $L$ used in this training stage, and it becomes:

[TABLE]

where $0<\alpha<1$, and it is a loss weighting coefficient (we used $\alpha=0.1$). This stage continues until the convergence or the termination condition is met (we restrict it to 25 epochs).

III Experiments

The objectives of our experiments are multi-fold. We verify the abilities of our unsupervised classification technique and compare it with other state-of-the-art clustering methods: $k$-means, where $k$ equals the number of target classes in the benchmark sets, and Gaussian mixture modeling (being a generalization of $k$-means which incorporates information about the covariance structure of the data and the centers of latent Gaussians), applied over the original (full) and reduced HSI. In the latter case, we reduce the dimensionality of the input HSI to match the size of our latent vector using (i) principal component analysis (PCA), (ii) independent component analysis (ICA), (iii) our sliding-window algorithm for simulating multispectral imagery from its hyperspectral counterpart, in which we generate the averaged band within a non-overlapping sliding window (S-MSI) [20], and (iv) our CAE. In this letter, the latent vector is of size $25$ (Fig 1). Also, we apply our CAE over reduced HSI and check the impact of the dimensionality reduction on its abilities (in this case, CAEs do not perform dimensionality reduction, and the latent vector is of the same size as the input vector). Finally, we compare unsupervised segmentation with our 1D-CNN [11] (Fig. 2) trained in a supervised manner over original and reduced benchmarks. Our study was divided into two experiments, over the available benchmarks (Section III-A), and a real-life HSI for which the ground-truth segmentation does not exist (Section III-B).

We exploited three most popular HSI benchmarks from the literature (http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes): (i) Salinas Valley (Sa), USA ($217\times 512$ pixels, AVIRIS sensor) showing different sorts of vegetation (16 classes, 224 bands, 3.7 m spatial resolution); (ii) Indian Pines (IP), USA ($145\times 145$, AVIRIS)—agriculture and forest (16 classes, 200 bands, 20 m); (iii) Pavia University (PU), Italy ($340\times 610$, ROSIS)—urban scenery (9 classes, 103 channels, 1.3 m). We also utilize the aerial hyperspectral observations acquired using the HyMap airborne sensor ($7982\times 512$, HyVista Corp. Pty Ltd., Australia, 126 bands with a wavelength resolution of 20 nm, 4.2 m) on Oct. 29, 2009. The study area was located in Mullewa, Western Australia (480 km2), and it is mainly used for the wheat, canola and lupin production. Although there were 30 test fields in which in-situ measurements had been performed (captured 1 m above a head of the wheat, eight measurements at each point; the measurement points are rendered in violet in Fig. 3), such data is not suitable for verifying segmentation algorithms, as we know the class label of an extremely small subset of all pixels. Hence, for Mullewa, we focus on qualitative analysis.

We use two clustering-quality measures to quantify the performance of the unsupervised techniques: normalized mutual information (NMI) and adjusted rand score (ARS). NMI is:

[TABLE]

where ${\rm MI}(A,B)=H(A)+H(B)-H(A,B)$ is the mutual information index quantifying the value of information shared between two random variables $A$ and $B$, $H(\cdot)$ denotes entropy, and $H(A,B)$ is the joint entropy between clusterings. ARS is:

[TABLE]

where $n$ is the number of objects (pixels) subjected to clustering, and $a$, $b$, $c$, and $d$ denote the number of data points placed in: the same group (cluster) in $A$ and $B$ ($a$), the same groups in $A$ and in different groups in $B$ ($b$), the same groups in $B$ and in different groups in $A$ ($c$), and in different groups in $A$ and in different groups in $B$ ($d$) [21]. Both NMI and ARS range from 0 to 1, where 1 means perfect score. Additionally, for 1D-CNN trained in a supervised manner using our balanced division into the training and validation sets (as presented in [11]; the B division), we report the average accuracy (AA), overall accuracy (OA), and the kappa scores $\kappa=1-\frac{1-p_{o}}{1-p_{e}}$, where $p_{o}$ and $p_{e}$ are the relative observed agreement, and hypothetical probability of chance agreement, respectively, and $-1\leq\kappa\leq 1$ ($\kappa=1$ is the perfect score). These scores were obtained over the entire input HSI to make them comparable with the unsupervised segmentation performed over the entire HSI (in both cases, however, we exclude the background pixels for which the class label is unknown), and averaged across 30 runs. Since the test set for 1D-CNN includes the training and validation examples, the results can be considered over-optimistic [11]. The deep networks were coded in Python 3.6, and the supervised training of 1D-CNN (ADAM, learning rate of $10^{-4}$, $\beta_{1}=0.9$, $\beta_{2}=0.999$) terminated, if after 25 epochs OA over the validation set (random subset of the training set) does not change. The experiments ran on NVIDIA GeForce RTX 2080.

III-A Experiment 1: Benchmark data

In this experiment, we compare 3D-CAE with other techniques over three HSI benchmarks—each unsupervised approach was executed exactly once, in order to understand its real-life applicability, where running algorithms multiple times is infeasible. Also, we performed Monte-Carlo cross-validation (repeated $30\times$) with balanced training and validation sets [11], and analyze the average supervised measures (AA, OA, and $\kappa$) obtained using 1D-CNN. For the sake of completeness, we report the unsupervised segmentation measures (NMI and ARS) for 1D-CNN as well—the entire scene (without background) was segmented. Since the test set includes both training and validation sets in this case (there is an “training-test information leak”), NMI and ARS may be considered over-optimistic for 1D-CNN.

In Table I, we gather the experimental results obtained over all sets. They show that 3D-CAE consistently delivers high-quality segmentation in all settings, with and without HSI reduction (in all cases, we decrease the feature dimensionality to $25$ to match the number of 3D-CAE embedded features). On the other hand, the dimensionality reduction is beneficial in the unsupervised setting, and leads to better clustering. It indicates that only a small portion of the entire spectrum conveys useful information about the captured materials within those HSI—exploiting the full spectrum makes segmentation much harder due to the curse of dimensionality (the best results were obtained using our S-MSI; Wilcoxon test, $p<0.001$). These observations are confirmed in Table II, where we report the ranking of all methods averaged across the benchmarks.

The execution time of all unsupervised techniques is reported in Table II. These times reflect only segmentation, without feature extraction for the methods run over reduced HSI (PCA took 1.17 s, 0.35 s, and 1.95 s for Sa, IP, and PU, respectively, ICA: 25.04 s, 1.28 s, and 56.68 s, S-MSI: 0.28 s, 0.04 s, 0.47 s, and 3D-CAE: 2843.27 s, 273.51 s, and 2975.91 s). Although 3D-CAE was significantly slower than other algorithms333The execution time of 3D-CAE over reduced HSI was very consistent with other deep learning-powered HSI segmentation techniques [17]., it retrieved consistently better segmentation (Table I). Also, we did not exploit early stopping for the clustering phase of 3D-CAE (it ran always for 25 epochs). This part of the training could have been terminated much earlier (as the training converged), which could have greatly reduced its execution time. It however requires further investigation.

III-B Experiment 2: Real-life data

In this experiment, we ran all unsupervised methods over a real-life hyperspectral scene. Since there is no ground-truth segmentation of the Mullewa dataset, we qualitatively compare the selected methods in Fig. 4. Here, we present the segmentations obtained using all unsupervised techniques over (i) full HSI, and (ii) reduced HSI (this reduction was performed with the approach which was the best over all benchmarks for the corresponding segmentation algorithm). We can appreciate that $k$-means and 3D-CAE give much more detailed segmentation (see example regions annotated with the white and black arrows in the GM visualization). It indicates that those regions are “heterogeneous” and manifest subtle spectral variations. This observation can trigger more detailed in-situ measurements (performed in precise locations), hence allow us to better understand the scanned regions and their critical characteristics. As previously, the execution time of 3D-CAE was much longer than other methods (Table II)—this issue can be tackled by more aggressive pre-processing (e.g., band selection), parallel GPU training or by applying early stopping conditions to both training phases of 3D-CAE.

IV Conclusion

We proposed a new deep learning-powered unsupervised HSI segmentation algorithm which exploits 3D convolutional autoencoders to learn embedded featues, and a clustering layer to segment an input image using the learned representation. Our experimental study, performed over benchmark and real-life HSI revealed that our approach delivers consistent and high-quality segmentation without any prior class labels. Such unsupervised techniques offer new possibilities to understand the acquired HSI—they can be used to: (i) enable practitioners to generate ground-truth HSI data in affordable time even for very large scenes (unsupervised segmentation of an input HSI would be reviewed and fine-tuned if necessary), (ii) perform anomaly detection within a captured region by analyzing unexpected heterogeneous parts of the segmentation map (e.g., a wheat farmland should be moderately homogeneous, and any deviation may be alarming), and to (iii) see beyond the current ground-truth HSI (Fig. 5). Although our method is computationally expensive, its execution time can be greatly decreased by the initial HSI reduction, applying early stopping conditions in both training phases, performing the parallel training (using multiple GPUs) and optimizing the hyper-parameters of the deep network architecture (e.g., decreasing the number of kernels)—it constitutes our current work.