# Feature Extraction for Hyperspectral Imagery: The Evolution From Shallow to Deep: Overview and Toolbox

Behnood Rasti, Danfeng Hong, Renlong Hang, Pedram Ghamisi, Xudong Kang, Jocelyn Chanussot, Jon Atli Benediktsson

<figure><img src="image_1.png" /></figure>

## To cite this version:

Behnood Rasti, Danfeng Hong, Renlong Hang, Pedram Ghamisi, Xudong Kang, et al. Feature Extraction for Hyperspectral Imagery: The Evolution From Shallow to Deep: Overview and Toolbox. IEEE geoscience and remote sensing magazine, 2020, 8 (4), pp.60-88. 10.1109/MGRS.2020.2979764. hal-03142175

HAL Id: hal-03142175

https://hal.science/hal-03142175v1

Submitted on 26 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Feature Extraction for Hyperspectral Imagery:
The Evolution from Shallow to Deep (Overview
and Toolbox)

Behnood Rasti, *Senior Member, IEEE*, Danfeng Hong, *Member, IEEE*, Renlong Hang, *Member, IEEE*, Pedram Ghamisi, *Senior Member, IEEE*, Xudong Kang, *Senior Member, IEEE*, Jocelyn Chanussot, *Fellow, IEEE*, and Jon Atli Benediktsson, *Fellow, IEEE*

**Abstract**—The final version of the paper can be found in IEEE Geoscience and Remote Sensing Magazine. Hyperspectral images provide detailed spectral information through hundreds of (narrow) spectral channels (also known as dimensionality or bands) with continuous spectral information that can accurately classify diverse materials of interest. The increased dimensionality of such data makes it possible to significantly improve data information content but provides a challenge to the conventional techniques (the so-called curse of dimensionality) for accurate analysis of hyperspectral images. Feature extraction, as a vibrant field of research in the hyperspectral community, evolved through decades of research to address this issue and extract informative features suitable for data representation and classification. The advances in feature extraction have been inspired by two fields of research, including the popularization of image and signal processing as well as machine (deep) learning, leading to two types of feature extraction approaches named shallow and deep techniques. This article outlines the advances in feature extraction approaches for hyperspectral imagery by providing a technical overview of the state-of-the-art techniques, providing useful entry points for researchers at different levels, including students, researchers, and senior researchers, willing to explore novel investigations on this challenging topic. In more detail, this paper provides a bird's eye view over shallow (both supervised and unsupervised) and deep feature extraction approaches specifically dedicated to the topic of hyperspectral feature extraction and its application on hyperspectral image classification. Additionally, this paper compares 15 advanced techniques with an emphasis on their methodological foundations in terms

This work is partially supported by Alexander von Humboldt (AvH) research grant.

B. Rasti is with the Machine Learning Group, Exploration Division, Helmholtz Institute Freiberg for Resource Technology, Helmholtz-Zentrum Dresden-Rossendorf, 09599 Freiberg, Germany (e-mail: behnood.rasti@gmail.com).

D. Hong is with Univ. Grenoble Alpes, CNRS, Grenoble INP, GIPSA-lab, 38000 Grenoble, France. (e-mail: hongdanfeng1989@gmail.com)

R. Hang is with the Jiangsu Key Laboratory of Big Data Analysis Technology, School of Automation, Nanjing University of Information Science and Technology, Nanjing 210044, China (e-mail:renlong_hang@163.com)

P. Ghamisi is with the Machine Learning Group, Exploration Division, Helmholtz Institute Freiberg for Resource Technology, Helmholtz-Zentrum Dresden-Rossendorf, 09599 Freiberg, Germany. (e-mail: p.ghamisi@gmail.com)

X. Kang is with the College of Electrical and Information Engineering, Hunan University, Changsha 410082, China, and also with the Key Laboratory of Visual Perception and Artificial Intelligence of Hunan Province, Changsha 410082, China (e-mail: xudong_kang@163.com).

J. Chanussot is with the Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, F-38000 Grenoble, France, also with the Faculty of Electrical and Computer Engineering, University of Iceland, Reykjavik 101, Iceland. (e-mail: jocelyn@hi.is)

J. A. Benediktsson is with the Faculty of Electrical and Computer Engineering, University of Iceland, 107 Reykjavik, Iceland (e-mail: benedikt@hi.is).

of classification accuracies. Furthermore, to push this vibrant field of research forward an impressive amount of codes and libraries is shared at: https://github.com/BehnoodRasti/HyFTech-Hyperspectral-Shallow-Deep-Feature-Extraction-Toolbox.

**Index Terms**—Classification, deep feature extraction, deep learning, dimensionality reduction, feature extraction, hyperspectral image, machine learning, shallow feature extraction, signal processing.

## I. INTRODUCTION

HYPERSPECTRAL imaging technology provides detailed spectral information by sampling the reflective portion of the electromagnetic spectrum covering a wide range from the visible region (0.4-0.7 µm) to the short-wave infrared (SWIR) region (almost 2.4 µm). Hyperspectral sensors can also characterize the emissive properties of objects by acquiring data in the range of the mid-wave and long-wave infrared regions, in hundreds of narrow contiguous spectral channels.

Detailed spectral information provided by hyperspectral sensors present both challenges and opportunities. For instance, hyperspectral images can be used to differentiate between different classes of interest with slightly different spectral characteristics [1]. However, most of the commonly used methods utilized for the analysis of gray scale, color, or multispectral images cannot be extended to analyse hyperspectral images due to several reasons, which will be detailed later in Section I.A.

The limited availability of training samples (which is a common issue in remote sensing) dramatically impacts the performances of supervised classification approaches due to the high dimensionality of hyperspectral images, which poses a problem for designing robust statistical estimations. Feature extraction can be used as a solution to address the aforementioned issue, which can be described as finding a set of vectors that represents an observation while reducing the dimensionality by transforming the input data linearly or nonlinearly to another domain and extract informative features in the new domain. The use of feature extraction techniques can be advantageous for a number of reasons which is illustrated in Fig. 1 and described in the following subsections.
<figure><img src="image_3.png" /><figcaption>Fig. 1. Feature extraction and its advantages for hyperspectral image analysis.</figcaption></figure>

A. Unique Geometrical and Statistical Properties of High Dimensional Data and the Need for Feature Extraction

Several studies (e.g., [2]–[4]) have demonstrated the unique geometrical, statistical, and asymptotical properties of high-dimensional data compared to RGB and multispectral images. These properties, which have been shown through some experimental and theoretical examples, clearly justify the reasons why most of the analytical approaches that are developed for RGB and multispectral images are not applicable to hyperspectral images [5]. Among those experimental examples, we can recall that (1) as dimensionality increases, the volume of a hypercube concentrates in corners, or (2) as dimensionality increases, the volume of a hypersphere concentrates in an outside shell. With respect to these examples, the following conclusions have been drawn:

* A high-dimensional feature space is almost empty, which indicates that multivariate data in $\mathbb{R}^p$ ($p$ represents the number of bands, spectral channels, or dimensions) can usually be represented in a lower dimensional space (referred to as subspace) without losing considerable information in terms of class separability [5].

* Since the high-dimensional feature space is almost empty (i.e., Gaussian distributed data have a tendency to concentrate in the tails while uniformly distributed data have a tendency to be concentrated in the corners), the density estimation of hyperspectral data for both Gaussian and uniform distributions become extremely challenging.

Fukunaga [6] claimed that there is a relation between the type of the classifier, the required number of training samples, and the number of input dimensions. As reported in [6], the required number of training samples is linearly related to the dimensionality for linear classifiers and to the square of the dimensionality for quadratic classifiers (e.g., the Gaussian maximum likelihood classifier [6]), and for nonparametric classifiers, the number of required samples exponentially increases as the dimensionality increases. Landgrebe showed a groundbreaking fact that too many spectral bands might have negative impacts in terms of expected classification performance [7]. When dimensionality increases, with a constant and limited number of training samples, a higher amount of statistics must be estimated. Thus, the accuracy of the statistical estimation

decreases although higher spectral dimensions increase the separability between the classes. This leads to a decrease in classification accuracies beyond an unknown number of bands. These problems are related to the *curse of dimensionality*, also known as *Hughes phenomenon* [8]. This finding was against the general understanding of hyperspectral data where it was wrongly believed that full dimensionality is always better than subspace in terms of classification accuracies.

The unique characteristics of high-dimensional data, as discussed above, have a pronounced impact on the performances of supervised classifiers [9], as they demand an adequate number of training samples, which is almost impossible to obtain in practice since the collection of such training samples is either expensive or time demanding. To address this issue, feature extraction-based dimensionality reduction is found to be effective.

B. Storage Systems and Processing Times

We are now in the era of massive data acquisitions. Statistics demonstrate that the cumulative volume of existing big data has been tremendously increased from 4.4 ZB to 44 ZB from 2013 to 2020¹. The EO community has also faced a similar trend because of the enormous volume and variety of data being generated by EO missions. For example, EnMAP (Environmental Monitoring and Analysis Program), which is a hyperspectral satellite mission, is planning to capture hyperspectral data with a maximum ground coverage of 5000 km × 30 km per day and the target revisit time of 4 days (±30°) with 512 Gbit onboard mass memory [10].

Feature extraction-based dimensionality reduction helps in data compression, which leads to the reduced storage space, faster transmission time, removing redundant features, reducing the storage space required, and fasten the required time for performing the same computations.

C. An Ever-Growing Relation between Machine learning and Feature Extraction

Fig. 2 illustrates the basic idea of a machine learning approach, which consists of feature extraction and classification.

<footnote n="1">¹https://www.emc.com/leadership/digital-universe/2014iview/executive-summary.htm</footnote>
In machine learning, users are requested to provide some guidelines for the machine (algorithm). These guidelines are usually provided by applying hand-crafted feature extraction approaches to provide informative features for the subsequent classifier. At the very beginning, each image pixel is regarded as a pattern, and its spectrum (i.e., a vector of different values of a pixel in different spectral channels) is considered as the initial set of features. This set of features, which are also known as spectral features, suffers from two important downsides: (1) These features are often redundant and (2) they do not consider the spatial dependencies of the adjacent pixels.

To address the first issue, a feature reduction (through feature extraction or selection) step can be applied, aiming at reducing the dimensionality of the input data (from $p_1$ dimensions in the original data to $p_2$ dimensions in a new feature space $p_2 < p_1$). This step can also be known as spectral feature extraction, which tries to preserve the key spectral information of the data by reducing the dimensionality and maximizing separability between classes. It is interesting to note that the second issue can also be addressed using feature extraction approaches. Please note that here the feature extraction step, which is also known as spatial feature extraction, is not aiming at reducing the dimensionality and instead is trying to model (extract) spatial contextual information suitable for the subsequent classification or object detection step and usually leads to an increase in the number of features. The simultaneous use of spectral and spatial features for hyperspectral data classification were studied in numerous works such as [5], [11]–[13].

Deep learning (as shown in Fig. 2), which is regarded as a subset of machine learning, tries to automatize the building blocks of the machine learning approaches (i.e., feature extraction and classification) by developing an end-to-end framework, which takes the input, performs automatic feature extraction and classification by considering the unique nature of the input data (instead of those hand-crafted feature extraction designs in machine learning), and outputs classification maps. It turns out that if an adequate amount of training data is supplied, deep learning approaches can outperform any other shallow machine learning approaches in terms of accuracies. Here, a question arises: Due to the fact that the number of available training data is often limited in the remote sensing community, would advanced deep learning-based approaches outperform their shallow alternatives in terms of accuracies? This question will be addressed in this paper.

Based on the above-mentioned descriptions, feature extraction is a key step in both machine learning and deep learning, whose concept has been evolved significantly through time from unsupervised to (semi-)supervised, from spectral or spatial to spectral-and-spatial, from manual to automatic, from hand-crafted to end-to-end, and from shallow to deep.

## D. Contributions

This article provides a detailed and organized overview of hyperspectral feature extraction techniques, categorized into two general sections: shallow feature extraction techniques (further categorized into supervised and unsupervised) and

<figure><img src="image_4.png" /><figcaption>Fig. 2. Feature extraction via machine learning and deep learning.</figcaption></figure>

deep feature extraction techniques. Each section provides a critical overview of the state of the art that is mainly rooted in the signal and image processing, statistical inference, and machine (deep) learning fields. Then, a few representative and advanced feature extraction approaches are chosen from each of the above-mentioned categories for further analysis and comparisons (mostly in terms of usefulness for classification). This article will, therefore, contribute to answering the following questions: When it comes to hyperspectral data in Earth observation, are deep learning-based feature extraction approaches better alternatives than their conventional (yet advanced) shallow feature extraction techniques? Which factors should be considered to design robust shallow and deep feature extraction techniques? In addition, to further promote this field of research, this paper is accompanied with a significant amount of codes and libraries for hyperspectral feature extraction, which is made publicly available at https://github.com/BehnoodRasti/HyFTech-Hyperspectral-Shallow-Deep-Feature-Extraction-Toolbox. Finally, several possible future directions are highlighted. To make the contribution of this paper clearer compared to the existing papers in the literature, here we provide a brief discussion. The work of [14] is dedicated to the evolution of discriminant analysis-based feature extraction models, which is a specific type of dimensionality reduction approaches. The work of [15] reviewed feature extraction and data mining works, which had been published mostly until 2012. Since 2012, however, many deep and shallow feature extraction approaches have been developed, which are critically reviewed and compared against each other in this work. The work of [16] focuses only on feature selection approaches while our proposed paper covers feature extraction techniques, and therefore, they complement each other.

## II. DATASETS AND NOTATIONS

### A. Datasets

1) *Indian Pines 2010:* This dataset (Fig. 3) is a very high resolution hyperspectral image (VHR HS) acquired by the ProSpecTIR system over near Purdue University, Indiana, between 24-25th of May 2010. In this paper, we use a subset
<figure><img src="image_5.png" /><figcaption>Fig. 3. Indian Pines 2010 datasets, from left to right, RGB composition, Training set, and Test set.</figcaption></figure>

of 445×750 pixels with 360 spectral bands. The dataset has the spatial resolution of 2 m and spectral width of 5 nm. The dataset contains 16 land cover classes shown in Fig. 3. The training and test sets used in this study have also been demonstrated in Fig. 3. Table I gives the number of samples, including training and test samples used in the experimental section.

<table><caption>TABLE I
INDIAN PINES 2010: THE NUMBER OF TRAINING SAMPLES, TEST SAMPLES, AND THE TOTAL NUMBER OF SAMPLES PER CLASS.</caption><thead><tr><th>Class No.</th><th>Class Name</th><th>Training Samples</th><th>Test Samples</th><th>Samples</th></tr></thead><tbody><tr><td>1</td><td>Corn-high</td><td>726</td><td>2661</td><td>3387</td></tr><tr><td>2</td><td>Corn-mid</td><td>465</td><td>1275</td><td>1740</td></tr><tr><td>3</td><td>Corn-low</td><td>66</td><td>290</td><td>356</td></tr><tr><td>4</td><td>Soy-bean-high</td><td>324</td><td>1041</td><td>1365</td></tr><tr><td>5</td><td>Soy-bean-mid</td><td>2548</td><td>35317</td><td>37865</td></tr><tr><td>6</td><td>Soy-bean-low</td><td>1428</td><td>27782</td><td>29210</td></tr><tr><td>7</td><td>Residues</td><td>368</td><td>5427</td><td>5795</td></tr><tr><td>8</td><td>Wheat</td><td>182</td><td>3205</td><td>3387</td></tr><tr><td>9</td><td>Hay</td><td>1938</td><td>48107</td><td>50045</td></tr><tr><td>10</td><td>Grass/Pasture</td><td>496</td><td>5048</td><td>5544</td></tr><tr><td>11</td><td>Cover crop 1</td><td>400</td><td>2346</td><td>2746</td></tr><tr><td>12</td><td>Cover crop 2</td><td>176</td><td>1988</td><td>2164</td></tr><tr><td>13</td><td>Woodlands</td><td>1640</td><td>46919</td><td>48559</td></tr><tr><td>14</td><td>Highway</td><td>105</td><td>4758</td><td>4863</td></tr><tr><td>15</td><td>Local road</td><td>52</td><td>450</td><td>502</td></tr><tr><td>16</td><td>Buildings</td><td>40</td><td>506</td><td>546</td></tr><tr><td></td><td>Total</td><td>10954</td><td>187120</td><td>198074</td></tr></tbody></table>

2) Houston University 2013: This dataset was acquired on June 23, 2012, by the Compact Airborne Spectrographic Imager (CASI) over the campus of the University of Houston and the neighboring urban area. The average height of the sensor was 5500ft. The data contain 349 × 1905 pixels with the spatial resolution of 2.5 m and 144 spectral bands ranging 0.38-1.05 µm. The dataset includes 15 classes of interests shown in Fig. 4. A color composite representation of the data and the corresponding training and test samples used in this study are shown in Fig. 4. The number of training

<figure><img src="image_6.png" /><figcaption>Fig. 4. Houston University 2013 datasets, from top to bottom, RGB composition, Training set, and Test set.</figcaption></figure>

and test samples for different classes of interests used in the experiments are given in Table II.

<table><caption>TABLE II
HOUSTON UNIVERSITY 2013: THE NUMBER OF TRAINING SAMPLES, TEST SAMPLES, AND THE TOTAL NUMBER OF SAMPLES PER CLASS.</caption><thead><tr><th>Class No.</th><th>Class Name</th><th>Training Samples</th><th>Test Samples</th><th>Samples</th></tr></thead><tbody><tr><td>1</td><td>Grass Healthy</td><td>198</td><td>1053</td><td>1251</td></tr><tr><td>2</td><td>Grass Stressed</td><td>190</td><td>1064</td><td>1254</td></tr><tr><td>3</td><td>Grass Synthetic</td><td>192</td><td>505</td><td>697</td></tr><tr><td>4</td><td>Tree</td><td>188</td><td>1056</td><td>1244</td></tr><tr><td>5</td><td>Soil</td><td>186</td><td>1056</td><td>1242</td></tr><tr><td>6</td><td>Water</td><td>182</td><td>143</td><td>325</td></tr><tr><td>7</td><td>Residential</td><td>196</td><td>1072</td><td>1268</td></tr><tr><td>8</td><td>Commercial</td><td>191</td><td>1053</td><td>1244</td></tr><tr><td>9</td><td>Road</td><td>193</td><td>1059</td><td>1252</td></tr><tr><td>10</td><td>Highway</td><td>191</td><td>1036</td><td>1227</td></tr><tr><td>11</td><td>Railway</td><td>181</td><td>1054</td><td>1235</td></tr><tr><td>12</td><td>Parking Lot 1</td><td>192</td><td>1041</td><td>1233</td></tr><tr><td>13</td><td>Parking Lot 2</td><td>184</td><td>285</td><td>469</td></tr><tr><td>14</td><td>Tennis Court</td><td>181</td><td>247</td><td>428</td></tr><tr><td>15</td><td>Running Track</td><td>187</td><td>473</td><td>660</td></tr><tr><td></td><td>Total</td><td>2832</td><td>12197</td><td>15029</td></tr></tbody></table>

3) Houston University 2018: This dataset was acquired on Feb. 16, 2017 by the hyperspectral imager CASI 1500 over the area of the University of Houston. The data covers the spectral range 380-1050 nm with 48 bands with the ground sampling distance of 1 m. In this article, we utilized the training portion of the whole data set, which was distributed by the Image Analysis and Data Fusion Technical Committee of the IEEE Geoscience and Remote Sensing Society (GRSS) and the University of Houston for the 2018 data fusion contest. It contains 601 × 2384 pixels and 20 land cover classes of interest shown in Fig. 5. The VHR RGB image is downsampled and shown in Fig. 5 together with the corresponding training and test samples used in this study. The number of training and test samples for different classes of interest used in the experiments are given in Table III.
<figure><img src="image_7.png" /></figure>

<table><caption>Fig. 5. Houston University 2018 datasets, from top to bottom, The VHR RGB Image (downsampled), Training set, and Test set.</caption><thead><tr><th></th><th>Healthy grass</th><th>Stressed grass</th><th>Synthetic grass</th><th>Evergreen Trees</th><th>Deciduous Trees</th></tr></thead><tbody><tr><th>Background</th><td></td><td></td><td></td><td></td><td></td></tr><tr><th>Soil</th><td>Water</td><td>Residential</td><td>Commercial</td><td>Road</td><td>Sidewalk</td></tr><tr><th>Crosswalk</th><td>Major Thoroughfares</td><td>Highway</td><td>Railway</td><td>Paved Parking Lot</td><td>Gravel Parking Lot</td></tr><tr><th>Cars</th><td>Trains</td><td>Seats</td><td></td><td></td><td></td></tr></tbody></table>

<caption>TABLE III</caption>

<table><caption>HOUSTON UNIVERSITY 2018: THE NUMBER OF TRAINING SAMPLES, TEST SAMPLES, AND THE TOTAL NUMBER OF SAMPLES PER CLASS.</caption><thead><tr><th>Class No.</th><th>Class Name</th><th>Training</th><th>Test</th><th>Sample</th></tr></thead><tbody><tr><td>1</td><td>Healthy grass</td><td>1458</td><td>8341</td><td>9799</td></tr><tr><td>2</td><td>Stressed grass</td><td>4316</td><td>28186</td><td>32502</td></tr><tr><td>3</td><td>Synthetic grass</td><td>331</td><td>353</td><td>684</td></tr><tr><td>4</td><td>Evergreen Trees</td><td>2005</td><td>11583</td><td>13588</td></tr><tr><td>5</td><td>Deciduous Trees</td><td>676</td><td>4372</td><td>5048</td></tr><tr><td>6</td><td>Soil</td><td>1757</td><td>2759</td><td>4516</td></tr><tr><td>7</td><td>Water</td><td>147</td><td>119</td><td>266</td></tr><tr><td>8</td><td>Residential</td><td>3809</td><td>35953</td><td>39762</td></tr><tr><td>9</td><td>Commercial</td><td>2789</td><td>220895</td><td>223684</td></tr><tr><td>10</td><td>Road</td><td>3188</td><td>42622</td><td>45810</td></tr><tr><td>11</td><td>Sidewalk</td><td>2699</td><td>31303</td><td>34002</td></tr><tr><td>12</td><td>Crosswalk</td><td>225</td><td>1291</td><td>1516</td></tr><tr><td>13</td><td>Major Thoroughfares</td><td>5193</td><td>41165</td><td>46358</td></tr><tr><td>14</td><td>Highway</td><td>700</td><td>9149</td><td>9849</td></tr><tr><td>15</td><td>Railway</td><td>1224</td><td>5713</td><td>6937</td></tr><tr><td>16</td><td>Paved Parking Lot</td><td>1179</td><td>10296</td><td>11475</td></tr><tr><td>17</td><td>Gravel Parking Lot</td><td>127</td><td>22</td><td>149</td></tr><tr><td>18</td><td>Cars</td><td>848</td><td>5730</td><td>6578</td></tr><tr><td>19</td><td>Trains</td><td>493</td><td>4872</td><td>5365</td></tr><tr><td>20</td><td>Seats</td><td>1313</td><td>5511</td><td>6824</td></tr><tr><td></td><td>Total</td><td>34477</td><td>470235</td><td>504712</td></tr></tbody></table>

## B. Notations

The observed HSI is denoted by $\mathbf{X} \in \mathbb{R}^{p \times n}$ where $p$ and $n$ are the number of spectral bands and pixels in each band, respectively. $d$ indicates the dimension of the feature space (the subspace). $\mathbf{X}_m \in \mathbb{R}^{p \times m}$ where $m < n$ denotes the matrix which contains the training samples. $\mathbf{y}_m \in \mathbb{R}^{1 \times m}$ denotes the vector which contains the class labels where $y_i \in \{1, 2, \dots, k\}$ and $k$ denotes the number of classes. $\mathbf{I}$ is the identity matrix and $\hat{\mathbf{X}}$ is the estimate of matrix $\mathbf{X}$. The Frobenius norm is denoted by $\|\cdot\|_F$. $\text{tr}(\mathbf{X})$ denotes the trace of matrix $\mathbf{X}$. The definitions of the symbols used in the paper are given in Table IV.

# III. SHALLOW FEATURE EXTRACTION TECHNIQUES

## A. Unsupervised Feature Extraction Techniques

Unsupervised feature extraction (UFE) often refers to the FE techniques which do not incorporate the knowledge of the ground (ground reference or labeled samples) to extract features. UFE techniques often rely on intrinsic characteristic of the HSI data such as geometric, spatial or spectral information to extract the features. Arguably, the main advantage of UFE compared with the other FE techniques is the lack of need for the training samples, which is of great importance in the case of remote sensing datasets. In this paper, four major UFE groups widely-used for HSI analysis are studied, which are categorized in the following subsections. Fig. 6 illustrates the graphical abstracts of those groups. Before explaining the UFE techniques in more details, we briefly refer to three groups of widely-used FE techniques which could also be assumed as UFE, however, will not be studied in details in this paper due to their specific applications. The first group includes a range of approaches such as normalized differential vegetation index (NDVI) and normalized differential water index (NDWI) which often rely on the knowledge of the characteristics of the sensors. The second group includes unmixing techniques, which could also be assumed as UFE techniques. They often exploit optimization techniques to show the fractions of materials existing in pixels based on some assumptions on the spectral signatures of the materials. Therefore, the final features extracted represent different materials in the scene at the sub-pixel level [17]. The third group includes an impressive number of approaches based on mathematical morphology, which hierarchically extract spatial and contextual information from the input image and usually leads to a significant increase in the number of features [18].

1) *Conventional Data Projection/Transformation Techniques*: Numerous UFE techniques fall into this category. The conventional techniques categorized in this group are often designed to linearly project or transform the data, $\mathbf{X}$, in a lower dimensional feature space (also called subspace) exploiting different non-local intrinsic characteristics of the hyperspectral dataset. The transformation can be given by

$$ \mathbf{Z} = \mathbf{V}^T \mathbf{X} \quad (1) $$

where $\mathbf{Z}$ is the projected data in the lower dimensional space and $\mathbf{V}$ is the transformation matrix or the bases for the subspace. Arguably, principal component analysis (PCA) [19] can
<table><caption>TABLE IV
THE DIFFERENT SYMBOLS USED IN THIS PAPER AND THEIR DEFINITION</caption><thead><tr><th>Symbols</th><th>Definition</th></tr></thead><tbody><tr><td>x<sub>i</sub></td><td>the <i>i</i>th entry of the vector <b>x</b>.</td></tr><tr><td>X<sub>ij</sub></td><td>the (<i>i</i>, <i>j</i>)th entry of the matrix <b>X</b>.</td></tr><tr><td>x<sub>i</sub></td><td>the <i>i</i>th column of the matrix <b>X</b>.</td></tr><tr><td>x<sub>(j)</sub></td><td>the <i>j</i>th row of the matrix <b>X</b>.</td></tr><tr><td>||<b>x</b>||<sub>0</sub></td><td><i>l</i><sub>0</sub>-norm of the vector <b>x</b> i.e. the number of nonzero entries.</td></tr><tr><td>||<b>x</b>||<sub>1</sub></td><td><i>l</i><sub>1</sub>-norm of the vector <b>x</b>, obtained by &sum;<sub>i</sub> |x<sub>i</sub>|.</td></tr><tr><td>||<b>x</b>||<sub>2</sub></td><td><i>l</i><sub>2</sub>-norm of the vector <b>x</b>, obtained by &radic;(&sum;<sub>i</sub> x<sub>i</sub><sup>2</sup>).</td></tr><tr><td>||<b>X</b>||<sub>1</sub></td><td><i>l</i><sub>1</sub>-norm of the matrix <b>X</b>, obtained by &sum;<sub>i,j</sub> |<b>X</b><sub>ij</sub>|.</td></tr><tr><td>||<b>X</b>||<sub>F</sub></td><td>Frobenius-norm of the matrix <b>X</b>, obtained by &radic;(&sum;<sub>i,j</sub> <b>X</b><sub>ij</sub><sup>2</sup>).</td></tr><tr><td>X&#775;</td><td>the estimate of the matrix <b>X</b>.</td></tr><tr><td>tr(<b>X</b>)</td><td>the trace of the matrix <b>X</b>.</td></tr><tr><td>||<b>X</b>||<sub>TV</sub></td><td>total variation norm of the matrix <b>X</b> is obtained by &sum;<sub>i</sub> TV(<b>x</b><sub>(i)</sub>)</td></tr></tbody></table>

be considered as the most conventional UFE technique which has been widely used for hyperspectral analysis [20]. PCA captures the maximum variance of the signal by projecting the signal on the eigenvectors of the covariance matrix (C) using

$$ \max_{\mathbf{V}} \frac{\mathbf{V}^T \mathbf{C} \mathbf{V}}{\mathbf{V}^T \mathbf{V}}. \qquad (2) $$

A widely used HSI UFE technique is the maximum noise fraction (MNF) [21] or noise adjusted principal components (NAPC) [22] which seek a projection in which the signal to noise ratio (SNR) is maximized. MNF uses the following optimization

$$ \max_{\mathbf{V}} \frac{\mathbf{V}^T \mathbf{C} \mathbf{V}}{\mathbf{V}^T \mathbf{C}_n \mathbf{V}}, \qquad (3) $$

where $C_n$ is the noise covariance matrix. Another conventional technique is independent component analysis (ICA) [23]. ICA assumes a linear mixture model of the non-Gaussian independent source signals and the mixing matrix which both are simultaneously estimated and therefore ICA is referred to as blind source separation. ICA has been also widely used for hyperspectral image analysis [24].

To cope with the nonlinearity of the HSI data, the kernel (nonlinear) versions of the aforementioned techniques i.e., kernel MNF [25], kernel ICA (KICA) [26], and kernel PCA (KPCA) [27] have been also proposed. By using the kernel trick the data are projected into a feature space where the inner product are defined using a kernel function. KICA and KPCA have been used as UFE techniques for change detection and classification in [28] and [29], respectively. In [30], discrete wavelet transformation (DWT) has been used for hyperspectral feature extraction. Since, DWT does not reduce the dimension, in [30], linear discriminant analysis (LDA) has been exploited to reduce the dimension.

2) Band Clustering/Splitting and Merging-based Techniques: The top right sub-figure from Fig. 6 shows the basic steps of band clustering and merging-based feature extraction methods. As shown in this figure, the core idea behind this group of methods is to split the spectral bands into several groups in which the spectral bands have very high correlation. Hence, the proposed techniques often use similarity and dissimilarity criteria to split the spectral bands into several

non-overlapping groups. By selecting or fusing the bands of each group, some representative bands or features of different groups are obtained. Furthermore, followed by the merging step, some band filtering and processing operations can be also used to further improve the discrimination of the resulting features. This group of techniques is often computationally cheap, and thus has been widely used in real applications. On the other hand, spectral information is often neglected by the methods in this category.

For the band clustering and merging two algorithms are proposed in [31]. The first algorithm selects discriminative bases by considering all the classes simultaneously, however, the second algorithm selects the best bases for a pair of classes at a time. In [32], a hierarchical clustering algorithm was introduced to split and cluster the hyperspectral bands where the representative band for each cluster is selected based on both a mutual information criterion and a divergence-based criterion. Another band clustering technique was proposed in [33] where the splitting is done by minimizing a mutual information criterion applied on averaged bands iteratively. Iterative algorithms were proposed in [34], for both splitting and merging the bands. The splitting procedure is done using the Pearson correlation coefficient between adjacent bands. Then, the merging is applied by averaging over the split bands.

Besides splitting/clustering and merging hyperspectral bands, another operation is to further improve the feature discrimination by band filtering or processing. For example, a hyperspectral feature extraction using image fusion and recursive filtering was given in [35] where the adjacent bands are fused by averaging and then recursive filtering was used for extracting spatial information. In [36], the intrinsic image decomposition is applied for processing the merged bands, which can effectively remove information that is not related to the material of different objects. After that, multiple improved versions of intrinsic decomposition-based band processing methods are developed [37], [38]. In [39], a relative total variation-based structure extraction method is applied for band processing, so as to construct multi-scale structural features which are robust to image noise.

3) Low-Rank Reconstruction-based Techniques: Low-rank reconstruction-based feature extraction techniques proposed
by Rasti et al. [40]–[43] are based on finding an orthogonal subspace by minimizing a constrained cost function. They exploit low-rank models and reconstruction-based optimization frameworks to extract features. The optimization frameworks take into account the prior knowledge of the data using different types of penalties. Due to the noise assumption in the low-rank model used, this group of FE techniques is robust to noise. They are often computationally expensive compared to groups 1 and 2 due to the iterative algorithms used to solve the (non-convex) optimization problem.

Wavelet-based sparse reduced rank regression (WSRRR) [41] applies the sparsity prior on the wavelet coefficients considering that the projected data on wavelet bases are sparse. WSRRR uses the model

$$ \mathbf{X} = \mathbf{V}^T \mathbf{Q} \mathbf{D}_2 + \mathbf{N}, \quad (4) $$

where $\mathbf{D}_2$ represents 2D wavelet bases, $\mathbf{X}$ is the observed HSI, $\mathbf{V}$ contains the orthogonal subspace bases, and $\mathbf{N}$ is the noise and model error. WSRRR simultaneously estimates the low-rank projection matrix and the wavelet coefficients $\mathbf{W}$ which minimizes

$$ (\hat{\mathbf{V}}, \hat{\mathbf{Q}}) = \arg \min_{\mathbf{Q}, \mathbf{V}} \frac{1}{2} \|\mathbf{X} - \mathbf{V}^T \mathbf{Q} \mathbf{D}_2\|_F^2 + \sum_{j=1}^d \lambda_j \|\mathbf{q}_{(j)}^T\|_1 \quad \text{s.t.} \quad \mathbf{V}^T \mathbf{V} = \mathbf{I}. \quad (5) $$

Note that the extracted features are given by $\hat{\mathbf{F}} = \hat{\mathbf{Q}}\mathbf{D}_2$.

To capture the spatial (neighboring) information, orthogonal total variation component analysis (OTVCA) has been proposed in [42] where the HSI is modeled as

$$ \mathbf{X} = \mathbf{V}^T \mathbf{F} + \mathbf{N}, \quad (6) $$

where matrix $\mathbf{F}$ contains the unknown features. OTVCA assumes that the hyperspectral features are spatially piece-wise smooth and, therefore, exploits the total variation (TV) penalty and simultaneously estimates $\mathbf{F}$ and $\mathbf{V}$ using

$$ \hat{\mathbf{F}} = \arg \min_{\mathbf{F}} \frac{1}{2} \|\mathbf{X} - \mathbf{V}^T \mathbf{F}\|_F^2 + \lambda \sum_{j=1}^d \mathrm{TV}(\mathbf{f}_{(j)}^T), \quad (7) $$

where

$$ \mathrm{TV}(\mathbf{x}) = \left\| \sqrt{(\mathbf{D}_h(\mathbf{x}))^2 + (\mathbf{D}_v(\mathbf{x}))^2} \right\|_1 $$

and $\mathbf{D}_v$ and $\mathbf{D}_h$ are the matrix operators to calculate the first order vertical and horizontal differences, respectively, of a vectorized image. Recently, sparse and smooth low-rank analysis (SSLRA) was proposed in [43] which models the HSI based on a combination of sparse and smooth features

$$ \mathbf{X} = \mathbf{V}^T (\mathbf{F} + \mathbf{S}) + \mathbf{N}, \quad (8) $$

where $\mathbf{F}$ and $\mathbf{S}$ contain smooth and sparse features, respectively. SSLRA extracts simultaneously the sparse features, $\mathbf{S}$, and the smooth ones, $\mathbf{F}$, by taking into account both sparsity and TV penalties

$$ (\hat{\mathbf{F}}, \hat{\mathbf{S}}, \hat{\mathbf{V}}) = \frac{1}{2} \|\mathbf{X} - \mathbf{V}^T (\mathbf{F} + \mathbf{S})\|_F^2 + \lambda_1 \|\mathbf{F}\|_{\mathrm{TV}} + \lambda_2 \|\mathbf{S}\|_1 \quad \text{s.t.} \quad \mathbf{V}^T \mathbf{V} = \mathbf{I}. \quad (9) $$

4) Graph Embedding and/or Manifold Learning Techniques: Considering the nonlinear characteristic of HSIs, this group of FE techniques aims to capture the data manifold through local geometric structure of neighboring pixels in the feature space. Fig. 6 (bottom right) demonstrates the concept of manifold learning FE techniques applied on the Swiss roll dataset. The pink line in the left figure shows the Euclidean distance between two data points in 3D space. It is clear that this line is not an effective metric to measure the similarity of the two points selected in the Swiss roll dataset. On the other hand, after unfolding the dataset which can be represented in 2D space, in the right figure, the Euclidean distance between two data points shown by the pink line is a better representative of the similarity of the two point in the dataset. The FE techniques categorized in this group are indeed designed to capture such a manifold while representing the data in a lower dimensional feature space.

Graph embedding or manifold learning FE techniques often include three main steps, 1) neighborhood pixel selection 2) weight selection, and 3) embedding construction. Isometric mapping (ISOMAP) [44], [45] is a global geometric nonlinear feature extraction. ISOMAP searches for geodesic distances between data points. It includes three main steps; 1) Constructing a neighborhood graph of the data points. 2) Computing the shortest path distances between all data points in the neighborhood graph. 3) constructing the lower dimensional embedding vectors which preserves the path distances in the neighborhood graph. Locally linear embedding (LLE) [46], Laplacian eigenmaps (LE) [47], and locality preserving projection (LPP) [48] are also geometric nonlinear feature extraction based on graph embedding. LLE constructs the embedding graph in three steps; 1) Selecting the neighbors for data points using the K nearest neighbors. 2) Compute the weights $\mathbf{A}_{ij}$ that linearly reconstruct the data points using their neighbors by minimizing the following constrained least-squares

$$ \min_{\mathbf{A}} \sum_{i=1}^{n} \left\| \mathbf{x}_i - \sum_{j \in \phi_i} \mathbf{A}_{ij} \mathbf{x}_j \right\|_2^2 \quad \text{s.t.} \quad \sum_{j \in \phi_i} \mathbf{A}_{ij} = \mathbf{1}, \quad (10) $$

where $\phi_i(\mathbf{x}_i)$ contains the neighborhood pixels selected for $\mathbf{x}_i$. We should note that the constrained weights estimated from (10) for every data point are invariant to rotations, rescalings, and translations of that data point and its neighbors and therefore they characterize the intrinsic geometric properties of each neighborhood. 3) Constructing the lower dimensional embedding vectors $\mathbf{y}$ by minimizing

$$ \min_{\mathbf{z}} \sum_{i=1}^{n} \left\| \mathbf{z}_i - \sum_{j \in \phi_i} \mathbf{A}_{ij} \mathbf{z}_j \right\|_2^2 \quad \text{s.t.} \quad \sum_{i=1}^{n} \mathbf{z}_i = \mathbf{0}, \\ \frac{1}{n} \sum_{i=1}^{n} \mathbf{z}_i \mathbf{z}_i^T = \mathbf{I}. \quad (11) $$

We should note that the reconstruction weights $\mathbf{A}_{ij}$ are fixed in minimization (11) and therefore the intrinsic geometric properties of the data with dimension $p$ are invariant to such a transformation into a lower dimension $d$.
Transform-based method

<figure><img src="image_8.png" /></figure>

Low rank reconstruction-based method

<figure><img src="image_9.png" /></figure>

Band Clustering and Merging-based method

<figure><img src="image_10.png" /></figure>

Manifold learning-based method

<figure><img src="image_11.png" /><figcaption>Fig. 6. Four major categories of unsupervised feature extraction methods.</figcaption></figure>

In [49], a general frame-work for graph embedding is given by

$$ \min_{\mathbf{z}} \sum_{i,j \in \phi_i} \|\mathbf{z}_i - \mathbf{z}_j\|_2^2 \mathbf{W}_{ij} \quad \text{s.t.} \quad \mathbf{ZBZ}^{\mathrm{T}} = \mathbf{I}, \quad (12) $$

or equivalently

$$ \min_{\mathbf{Z}} \operatorname{tr}(\mathbf{Z}(\mathbf{D} - \mathbf{W})\mathbf{Z}^{\mathrm{T}}) = \min_{\mathbf{Z}} \operatorname{tr}(\mathbf{ZLZ}^{\mathrm{T}}) \quad \text{s.t.} \quad \mathbf{ZBZ}^{\mathrm{T}} = \mathbf{I}, \quad (13) $$

where $\mathbf{L} = \mathbf{D} - \mathbf{W}$ denotes the Laplacian matrix of the undirected weighted graph $G = \{\mathbf{X}, \mathbf{W}^p\}$ (where $\mathbf{X}$ is the vertex set and $\mathbf{W} \in \mathbb{R}^{n \times n}$ is the similarity matrix), $\mathbf{D}$ is a diagonal matrix where its entries are given by

$$ \mathbf{D}_{ii} = \sum_{j \neq i} \mathbf{W}_{ij}, \forall i. \quad (14) $$

The diagonal matrix $\mathbf{B}$ is for the scale normalization and might also be the Laplacian matrix of a penalty graph such as $G^p = \{\mathbf{X}, \mathbf{W}^p\}$. We should note that the vertices of $G^p$ and $G$ (i.e., $\mathbf{X}$) are the same while the similarity matrix ($\mathbf{W}^p$) corresponds to the similarity characteristics suppressed in the lower dimensional feature space ($\mathbf{B} = \mathbf{L}^p = \mathbf{D}^p - \mathbf{W}^p$, see [49]). LLE can be reformulated using graph embedding mentioned above with similarity matrix $\mathbf{W}_{ij} = \mathbf{A}_{ij} + \mathbf{A}_{ij}^T - \mathbf{A}_{ij}^T \mathbf{A}_{ij}$ if $j \in \phi_i$ otherwise $\mathbf{W}_{ij} = 0$ and $\mathbf{B} = \mathbf{I}$ [49]. ISOMAP, LE, and LPP can be also formulated using graph embedding [49]. In the viewpoint of graph embedding the main differences between these FE techniques are the selection of the matrices $\mathbf{W}$ and $\mathbf{B}$. For instance, LE and LPP use the Gaussian function with the standard deviation $\sigma$ to choose the similarity matrix as

$$ \mathbf{W}_{ij} = \begin{cases} \exp \frac{-\|\mathbf{x}_i - \mathbf{x}_j\|_2^2}{2\sigma^2}, & \forall i, j \in \phi_i(\mathbf{x}_i) \\ 0, & \text{otherwise.} \end{cases} \quad (15) $$

We should note that the techniques categorized in this group are assumed as supervised FE methods when they are applied only on the training samples. This is common in the case of

HSI due to the large volume of the image which makes the algorithm computationally very expensive. In the following section we will discuss how the ground reference (training samples) can be used to construct the edge matrix, $\mathbf{W}$, and therefore those techniques are considered as SFE.

## B. Supervised Feature Extraction Techniques

Unlike unsupervised FE techniques that rely on modeling various prior assumptions of hyperspectral data supervised methods are capable of extracting class-separable features more effectively, owing to the use of label information. Over the past few decades, some seminal models have been widely developed and applied to perform supervised feature extraction on HSIs, which can be roughly categorized into two streams: *Subspace learning (SL)-based* and *band selection (BS)-based* approaches.

Different from the hand-crafted features [50], the SL-based approaches learn to extract the low-dimensional representation from the data by formulating different supervised rules in view of label information. There are some typical methods in SL, including linear discriminant analysis (LDA) [51], matrix discriminant analysis (MDA) [52], decision boundary FE (DBFE) [53], etc. While the latter, BS, which aims at screening out the representative and informative spectral bands, is unfolded with mutual information-based BS [54], rough set, and fuzzy C-means [55], to name a few. To further enhance the class separability, a large number of extended methods have been successfully proposed in recent years, which are subspace LDA (SLDA) [56], regularized LDA [57], local fisher's discriminant analysis (LFDA) [58], feature space discriminative analysis (FSDA) [59], rough-set-based BS [60], and FE with local spatial modeling [61].

According to the powerful learning ability of SL methods compared to that of BS-based strategies, we rather focus on reviewing the SL-related FE techniques, in which two main streams – discriminant analysis FE (DAFE) and regression-induced representation learning (RIRL) – are emphatically
<figure><img src="image_12.png" /><figcaption>Fig. 7. The illustration for supervised FE with four different categories. The obvious differences lie in the use form of label information and learning strategies, i.e., LDA: Fishers rule, CGDA: Affinity matrix, LSDR: Labels, JPlay: Joint use of affinity matrix and one-hot encoded labels.</figcaption></figure>

investigated and compared by clarifying their similarities and differences as well as pros and cons, as briefly illustrated in Fig. 7.

1) Discriminant Analysis Feature Extraction (DAFE): Generally speaking, DAFE seeks to find an optimal projection or transformation matrix $\mathbf{P} \in \mathbb{R}^{p \times d}$ ($d$ is the dimension of the to-be-estimated subspace) by optimizing certain class-relevant separation criterion associated with the label information. In this process, the estimated subspace $\mathbf{Z} \in \mathbb{R}^{d \times n}$ that consists of a series of vector $\mathbf{z}_i$ can be obtained by projecting the samples $\mathbf{X}_m = \{\mathbf{x}_i\}_{i=1}^m \in \mathbb{R}^{p \times m}$ onto a decision boundary, which can be generally expressed as $\mathbf{Z} = \mathbf{P}^T \mathbf{X}$. Each vector $\mathbf{z}_i$ in $\mathbf{Z}$ can be collected by $\mathbf{P}^T \mathbf{x}_i$. Depending on the different types of label embedding, DAFE can be subdivided into LDA and its variants, graph-based discriminant analysis (GDA) and its extensions, and kernelized discriminant analysis (KDA).

▷ LDA and Its Variants: The traditional LDA linearly transforms the original data into a discriminative subspace by maximizing the Fisher's ratio in the form of generalized Rayleigh quotient, that is, minimizing the intra-class scatter and maximizing inter-class scatter simultaneously. Given a pair-wise training set $\{(\mathbf{x}_1, \mathbf{y}_1), \dots, (\mathbf{x}_i, \mathbf{y}_i), \dots, (\mathbf{x}_m, \mathbf{y}_m)\}$, the objective function of multi-class LDA to estimate the linear projection matrix $\mathbf{P}$ can be written as follows:

$$ \max_{\mathbf{P}} \frac{\mathrm{tr}(\mathbf{P}^\mathrm{T}\mathbf{S}_b\mathbf{P})}{\mathrm{tr}(\mathbf{P}^\mathrm{T}\mathbf{S}_w\mathbf{P})}, \qquad (16) $$

where $\mathbf{S}_w$ and $\mathbf{S}_b$ are defined as the within-class scatter matrix and the between-class scatter matrix, respectively. With the constraint of $\mathbf{P}^\mathrm{T}\mathbf{S}_w\mathbf{P} = \mathbf{I}$, the optimization problem in (16) can be equivalently converted to one of $\mathbf{S}_b\mathbf{P} = \lambda\mathbf{S}_w\mathbf{P}$ by introducing the Lagrange multiplier $\lambda$. The close-form solution to the simplified optimization problem can be deduced by a generalized eigenvalues decomposition (GED).

Due to the sensitivity to complex high-dimensional noises caused by the environmental and instrumental factors and the availability of labeled samples, the original LDA inevitably suffers from an ill-posed statistical degradation, especially in the case of small-scale samples. The degraded reasons mainly lie in the singularity of the two scatter metrics ($\mathbf{S}_w$ and $\mathbf{S}_b$), thereby easily leading to the overfitting problem. To improve the stability and generalization, the regularized LDA was proposed by additionally adding a $l_2$-norm constraint on $\mathbf{S}_w$ parameterized by $\gamma$ as $\mathbf{S}_w^{\mathrm{reg}} = \mathbf{S}_w + \gamma\mathbf{I}$. By replacing $\mathbf{S}_w$ in (16) with the regularized $\mathbf{S}_w^{\mathrm{reg}}$, the solution in the regularized LDA can be still obtained by the GED solver.

Considering the local neighborhood relations between samples in the process of model learning, LFDA breaks through the bottleneck of those LDA-based methods by assuming that the data are distributed in the nonlinear manifolds rather than a homogeneous Gaussian space. For this purpose, LFDA is capable of effectively excavating the locally underlying structure of the data that lies in the real world. Essentially, LFDA can be regarded as a weighted LDA by locally weighing $\mathbf{S}_w$ and $\mathbf{S}_b$ matrices. Therefore, the two modified scatter matrices, denoted as $\tilde{\mathbf{S}}_w$ and $\tilde{\mathbf{S}}_b$, can be formulated as

$$ \begin{aligned} \tilde{\mathbf{S}}_w &= \frac{1}{2} \sum_{i=1}^{m} \sum_{j=1}^{m} \mathbf{W}_{ij}^w (\mathbf{x}_i - \mathbf{x}_j)(\mathbf{x}_i - \mathbf{x}_j)^\mathrm{T}, \\ \tilde{\mathbf{S}}_b &= \frac{1}{2} \sum_{i=1}^{m} \sum_{j=1}^{m} \mathbf{W}_{ij}^b (\mathbf{x}_i - \mathbf{x}_j)(\mathbf{x}_i - \mathbf{x}_j)^\mathrm{T}, \end{aligned} \qquad (17) $$

where the two weights ($\mathbf{W}^w$ and $\mathbf{W}^b$) denote the sample-wise similarities. There are several widely-used strategies in calculating such a similarity matrix symbolized by $\mathbf{W}$. A simple yet effective one is given by $\mathbf{W}_{ij} = 1$, if $\mathbf{x}_j \in \phi_k(\mathbf{x}_i)$, where $\phi_k(\mathbf{x}_i)$ represents the k-nearest-neighbor of $\mathbf{x}_i$; otherwise,
$W_{ij} = 0$. Another commonly-used one was constructed based on the radial basis function (RBF) with a standard derivation of $\sigma$, as defined in (15). Please refer to [62]–[64] that might be useful for those who are interested in more types of $W$.

Similar to SLDA that first projects the original data into a subspace and then LDA is performed in the transformed subspace, FSDA starts with maximizing the between-spectral scatter matrix ($S_f$) to enhance the differences along the spectral dimension, and similarly, the LDA is further used for extracting the representations of class separability from the feature domain. In the first step, let $\mu_{i,j}$ be the average value of the j-th class and the i-th spectral band, then we have the definition of $S_f$ as follows:

$$S_f = \frac{1}{2} \sum_{i=1}^{p} (\mathbf{h}_i - \bar{\mathbf{h}})(\mathbf{h}_i - \bar{\mathbf{h}})^{\mathrm{T}}, \quad (18)$$

where $\mathbf{h}_i = [\mu_{i,1}, \mu_{i,2}, ..., \mu_{i,k}]$ is the spectral representation in the feature space and $\bar{\mathbf{h}} = \frac{1}{p} \sum_{i=1}^{p} \mathbf{h}_i$. The primary transformation ($\mathbf{P}_f$) that aims at improving the spectral discriminant can be estimated with maximizing the trace term of $S_f$ as

$$\max_{\mathbf{P}_f} \mathrm{tr}(\mathbf{P}_f^{\mathrm{T}} \mathbf{S}_f \mathbf{P}_f). \quad (19)$$

Using the obtained $\mathbf{P}_f$, the latent representation in the feature space $\mathbf{g}_i = \mathbf{P}_f^{\mathrm{T}} \mathbf{h}_i$, $i = 1, 2, ..., p$ can be further fed into the next-step LDA.

**▷ GDA and Its Extensions:** Before revisiting the GDA methods, we first introduce and formulate the general graph embedding (GGE) framework presented in [49] with Eq. (12). Obviously, the extracted features $\mathbf{Z}$ in the GGE framework are determined by the construction of $\mathbf{W}$ to a great extent. Thus, we will highlight several kinds of representative affinity matrices corresponding to the different graph embedding approaches, i.e. LDA, LE [47] and its linearized LPP [48], LLE [46], sparse GDA (SGDA) [65], and collaborative GDA (CGDA) [66]. Fig. 8 visualizes the affinity matrices given by five different strategies in a four-class case selected from Houston 2013 dataset.

**LDA-like affinity matrix:** In essence, LDA is vested in a special case of GGE framework with $\mathbf{D}^{(\mathrm{LDA})} = \mathbf{I}$, whose affinity matrix can be represented as

$$\mathbf{W}_{ij}^{(\mathrm{LDA})} = \begin{cases} 1/N_k, & \text{if } \mathbf{x}_i \text{ and } \mathbf{x}_j \in C_k; \\ 0, & \text{otherwise,} \end{cases} \quad (20)$$

where $N_k$ is the number of samples belonging to k-th class.

**LPP or LE-based affinity matrix:** One is to be constructed in kernel space with a higher dimension via similarity measurement, i.e. extensively using (15).

**LLE-based affinity matrix:** Different from the hand-crafted graph, LLE reconstructs each given sample with its k-nearest neighbors by exploiting the linear regression techniques [67], [68]. As a result, the reconstruction coefficients ($\mathbf{A}$) can be obtained by solving the optimization problem of (10). With the known $\mathbf{A}$, it is straightforward to derive the needful affinity matrix, denoted as $\mathbf{W}^{(\mathrm{LLE})}$,

$$\mathbf{W}_{ij}^{(\mathrm{LLE})} = \begin{cases} \mathbf{A}_{ij} + \mathbf{A}_{ij}^{\mathrm{T}} - \mathbf{A}_{ij} \mathbf{A}_{ij}^{\mathrm{T}}, & \text{if } \mathbf{x}_j \in \phi_k(\mathbf{x}_i); \\ 0, & \text{otherwise,} \end{cases} \quad (21)$$

thereby inducing the Laplacian matrix as $\mathbf{L}^{(\mathrm{LLE})} = \mathbf{D}^{(\mathrm{LLE})} - \mathbf{W}^{(\mathrm{LLE})} = (\mathbf{I} - \mathbf{A})^{\mathrm{T}}(\mathbf{I} - \mathbf{A})$.

**SGDA and CGDA-guided affinity matrix:** Similarly to LLE, the affinity matrix can be estimated using the data-driven representation learning, i.e., sparse and collaborative representations [69]–[71]. Accordingly, the two learning strategies can be equivalent to respectively solving the constrained $l_1$-norm optimization problem:

$$\min_{\mathbf{W}} \|\mathbf{W}\|_1 \quad \text{s.t.} \quad \|\mathbf{X}_m \mathbf{W} - \mathbf{X}_m\|_F^2 \le \epsilon, \quad (22)$$

and the $l_2$-norm optimization problem:

$$\min_{\mathbf{W}} \|\mathbf{W}\|_F^2 \quad \text{s.t.} \quad \|\mathbf{X}_m \mathbf{W} - \mathbf{X}_m\|_F^2 \le \epsilon. \quad (23)$$

The aforementioned affinity matrices can be unified to the GGE framework of (12).

In addition to SGDA and CGDA (the two baselines), Huang et al. [72] learned a set of sparse coefficients on manifolds and then preserved the sparse manifold structure in the embedded space. The work in [73] extended the existing SGDA to the spatial-spectral graph embedding to address the issues of the spatial variability and spectral multimodality. With the embedding of the intrinsic geometric structure of the data, a Laplacian regularizer CGDA [74] was developed to further improve the graph's confidence. Li et al. [75] simultaneously integrated the sparsity and low-rankness into the graph for capturing a more robust structure of the data locally and globally. Furthermore, Pan et al. [76] further improved the above work by unfolding the HSI data with the form of a tensor.

**▷ KDA:** In reality, the HSI usually exhibits a highly nonlinear data distribution, which may result in difficulties to effectively identify the materials. The solution to this issue is making use of a so-called kernel trick [77] that can map the data of the input space into a new Hilbert space with a higher feature dimension. In the kernel-induced space, the complex nonlinearity of the HSI can be well analyzed in a linearized system. Comparatively, the input to KDA is an inner product of original data pairs, defined as $k(\mathbf{x}_i, \mathbf{x}_j)$ which can be given by (15). By introducing the kernel Gram matrix $\mathbf{K}$ with $\mathbf{K}_{i,j} = \mathbf{\Phi}(\mathbf{x}_i)^{\mathrm{T}} \mathbf{\Phi}(\mathbf{x}_j) = k(\mathbf{x}_i, \mathbf{x}_j)$, most of previous LDA-based methods can be simply extended to the corresponding kernelized versions, i.e. KLDA and KLFDA can calculate their projections $\mathbf{P}$ by solving a GED problem of

$$\mathrm{KLKP} = \lambda(\mathrm{KBK} + \gamma\mathbf{I})\mathbf{P}. \quad (24)$$

Note that $\mathbf{B} = \mathbf{I}$ in KLDA, while $\mathbf{L} = \mathbf{L}_w$ and $\mathbf{B} = \mathbf{L}_b$ are computed by $\mathbf{D}_w - \mathbf{W}_w$ and $\mathbf{D}_b - \mathbf{W}_b$ in the kernel space, respectively, for KLFDA. Furthermore, for KSGDA and KCGDA, the main difference lies in the computation of the adjacency matrix, which can be performed in the kernel space by solving the general kernel coding problem as follows:

$$\min_{\mathbf{W}} \Omega(\mathbf{W}) \quad \text{s.t.} \quad \|\mathbf{\Phi}(\mathbf{X}_m)\mathbf{W} - \mathbf{\Phi}(\mathbf{X}_m)\|_F^2 \le \epsilon, \quad (25)$$

where $\Omega(\mathbf{W})$ can be selected to be either sparsity-prompting term $\|\mathbf{W}\|_1$ of KSGDA or dense (or collaborative) term $\|\mathbf{W}\|_F^2$ of KCGDA. In [78] and [74], the solutions in (25) have been theoretically guaranteed in the same way by solving the
<figure><img src="image_13.png" /><figcaption>Fig. 8. A four-class showcase for affinity matrices (W) with respect to five different approaches, where the connectivity (or edge) of W is computed within each class.</figcaption></figure>

problems (22) and (23) using the alternating direction method of multiplier (ADMM) [79] and least-square regression with Tikhonov regularization [80], respectively.

## 2) Regression-induced Representation Learning (RIRL):

RIRL provides a new insight from the regression's point of view to model the FE behavior by bridging the training samples with the corresponding labels rather than indirectly using the label information in the form of graph or affinity matrix in DAFE-based methods.

### ▷ Least-Squares Dimension Reduction (LSDR):
We begin with sliced inverse regression (SIR) [81], which is a landmark in supervised FE techniques. It assumes that the pair-wise data $\{(x_i, y_i)\}_{i=1}^m$ are conditionally independent on the to-be-estimated subspace features $\{z_i\}_{i=1}^m$, formulated as $(X \perp Y) | Z$. Following this rule, the LSDR proposed by Suzuki and Sugiyama [82] attempts to find a maximizer of the squared-loss mutual information (SMI) to satisfy the previously mentioned independence assumption. The projections $\mathbf{P}$ for LSDR can be searched by optimizing the following maximization problem:

$$ \max_{\mathbf{P}} \text{SMI}(\mathbf{Z}, \mathbf{Y}) \quad \text{s.t.} \quad \mathbf{P}\mathbf{P}^\mathrm{T} = \mathbf{I}, \qquad (26) $$

and the SMI to measure a statistical dependence between two discrete variables is defined as

$$ \text{SMI}(\mathbf{Z}, \mathbf{Y}) = \sum_{\mathbf{z} \in \mathbf{Z}} \sum_{\mathbf{y} \in \mathbf{Y}} p(\mathbf{z})p(\mathbf{y}) \left( \frac{p(\mathbf{z}, \mathbf{y})}{p(\mathbf{z})p(\mathbf{y})} - 1 \right)^2, \quad (27) $$

where $p(\bullet)$ is the probability distribution function.

### ▷ Least-Squares Quadratic Mutual Information (LSQMI):
Limited by the sensitivity of MI to outliers, authors of [83] designed a more robust LSQMI with the basis of QMI criterion, hence let us define the QMI as

$$ \text{QMI}(\mathbf{Z}, \mathbf{Y}) = \sum_{\mathbf{z} \in \mathbf{Z}} \sum_{\mathbf{y} \in \mathbf{Y}} (p(\mathbf{z}, \mathbf{y}) - p(\mathbf{z})p(\mathbf{y}))^2. \quad (28) $$

Similarly, we solve (26)-like optimization problem by replacing SMI with QMI.

### ▷ Least-Squares QMI Derivative (LSQMID):
Due to the difficulty in accurately computing the derivative of QMI estimator, LSQMI was further extended to a computationally effective LSQMID by estimating the derivative of QMI instead of QMI itself [84]. In this work, authors have demonstrated a more accurate and efficient derivative computation of QMI.

### ▷ Joint & Progressive Learning Strategy (JPlay):
Another MI-free estimation group is latent subspace learning (LSL).

One representative LSL performs FE and classification simultaneously in joint learning (JL) fashion [85]. With an expected output $\Theta X_m$, the process can be modeled as

$$ \min_{\mathbf{P}_k, \Theta} \|\mathbf{Y}_l - \mathbf{P}_k \Theta \mathbf{X}_m\|_F^2 + \frac{\alpha}{2} \|\mathbf{P}_k\|_F^2 \quad \text{s.t.} \quad \Theta \Theta^\mathrm{T} = \mathbf{I}, \quad (29) $$

where $\mathbf{Y}_l \in \mathbb{R}^{k \times m}$ and $\Theta \in \mathbb{R}^{d \times p}$ are defined as the one-hot encoded label matrices and the latent subspace projections, respectively. $\mathbf{P}_k \in \mathbb{R}^{k \times d}$ denotes the regression matrix that connects the learned subspace and the label information. $\mathbf{Y}_l$ can be formulated as

$$ \mathbf{Y}_l = \begin{bmatrix} 1 & 0 & \dots & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & 0 & \dots & 1 \end{bmatrix} \begin{array}{c} 1 \\ 2 \\ \vdots \\ j \\ \vdots \\ k \end{array} \quad (30) $$

In [85], the model's solution has been proven to be a closed-form. Moreover, the work in [86] explored a LDA-like graph as a regularizer to learn a spectrally discriminative feature representation, thus (29) becomes

$$ \min_{\mathbf{P}_k, \Theta} \|\mathbf{Y}_l - \mathbf{P}_k \Theta \mathbf{X}_m\|_F^2 + \frac{\alpha}{2} \|\mathbf{P}_k\|_F^2 + \frac{\beta}{2} \mathrm{tr}(\Theta \mathbf{X}_m \mathbf{L} \mathbf{X}_m^\mathrm{T} \Theta^\mathrm{T}) \quad \text{s.t.} \quad \Theta \Theta^\mathrm{T} = \mathbf{I}. \quad (31) $$

Beyond the JL-based models, Hong et al. [87] established a novel multi-layered regression framework by following a joint and progressive learning strategy (JPlay). With the layer-wise auto-reconstruction mechanism effectively against spectral variabilities caused by complex noises and atmospheric effects, the linearized JPlay breaks through the performance bottleneck of traditional linear methods. More specifically, we have the resulting model in the following

$$ \begin{align} & \min_{\mathbf{P}_k, \{\Theta_l\}_{l=1}^q} \|\mathbf{Y}_l - \mathbf{P}_k \Theta_q \dots \Theta_1 \mathbf{X}_m\|_F^2 + \frac{\alpha}{2} \|\mathbf{P}_k\|_F^2 \nonumber \\ & \qquad + \frac{\beta}{2} \sum_{l=1}^q \mathrm{tr}(\Theta_l \mathbf{X}_{l-1} \mathbf{L} \mathbf{X}_{l-1}^\mathrm{T} \Theta_l^\mathrm{T}) \nonumber \\ & \qquad + \frac{\gamma}{2} \sum_{l=1}^q \|\mathbf{X}_{l-1} - \Theta_l^\mathrm{T} \Theta_l \mathbf{X}_{l-1}\|_F^2, \tag{32} \\ & \text{s.t.} \quad \mathbf{X}_l = \Theta_l \mathbf{X}_{l-1}, \quad \mathbf{X}_l \ge 0, \quad \|\mathbf{x}_i\|_2 \le 1, \nonumber \end{align} $$
where the soft constraint $||x_i||_2 \le 1$ can be used to relax the orthogonality. It is worth noting that such JL-based strategy can clearly tell the model which features are positive to classification task, owing to the joint strategy of feature extraction and classification.

## IV. DEEP FEATURE EXTRACTION TECHNIQUES

Shallow feature extraction techniques often require careful engineering and domain knowledge of experts, which limit their applications. Different from them, deep learning techniques aim at automatically learning high-level features from raw data in a hierarchical fashion. These features are more discriminative, abstract, and robust than shallow ones. Due to their powerful feature representation ability, deep learning techniques have been widely used to extract features from hyperspectral images in recent years [88], [89]. Among various deep learning models, autoencoders (AEs), convolutional neural networks (CNNs), and recurrent neural networks (RNNs), shown in Fig. 9, are the most popular ones. In this section, we will present these models and their applications to hyperspectral feature extraction.

### A. AEs

As demonstrated in Fig. 9, AE is mainly comprised of two modules: Encoder and decoder. Encoder maps the input vector $x$ into a hidden space $h$, while decoder aims at getting a reconstruction result $\hat{x}$ of the original input from $h$. These processes can be formulated as

$$ 
\begin{aligned}
\mathbf{h} &= f(\mathbf{W}_1\mathbf{x} + \mathbf{b}_1), \\
\hat{\mathbf{x}} &= f(\mathbf{W}_2\mathbf{h} + \mathbf{b}_2),
\end{aligned} 
\quad (33) 
$$

where $\mathbf{W}_1$ and $\mathbf{W}_2$ denote the weights connecting the input layer to the hidden layer and the hidden layer to the output layer, respectively, $\mathbf{b}_1$ and $\mathbf{b}_2$ represent the biases of the hidden units and the output units, respectively, and $f$ is a nonlinear activation function. The training of AE is to minimize the residual between $\mathbf{x}$ and $\hat{\mathbf{x}}$. Once trained, the decoder is deleted and the hidden layer $h$ is considered as a feature representation of $\mathbf{x}$. In order to extract deep features, several AEs are often stacked together, generating a stacked AE (SAE) model. For SAE, the hidden layer in the preceding AE will be used as the input of the subsequent AE.

SAE is perhaps the earliest deep model used to extract features of hyperspectral images [90]. One typical benefit of SAEs is that each AE inside the network can be pre-trained using both labeled and unlabeled samples, thus providing better initial values for network parameters compared to the random initialization. After layer-wise pre-training, the fine-tuning of only a few layers can acquire satisfactory discriminant features. This training method is capable of alleviating overfitting problem when there only exist small numbers of training samples in hyperspectral images. In [91], the spectral information for each pixel was considered as a vector, and fed into an SAE model to extract deep features. These features extracted by SAEs can also be generalized from one image to another image, which was validated in [92] and [93].

In order to extract spatial features of each pixel, one often needs to select a local patch or cube centered at the pixel, and then input it into a feature extraction model. Since the inputs of SAEs are vectors, it is difficult to directly process patches or cubes. In [91] and [92], the local cubes from the first principal components of hyperspectral images were firstly reshaped into vectors, and then fed into SAEs to extract spatial features. In [94] and [95], Gabor features and extended morphological attribute profiles (i.e., the joint use of shallow and deep feature extraction methods) were used as the inputs of SAEs, making the network easier to extract high-level spatial features. After the extraction of spectral and spatial features, they can be easily concatenated together to generate a spectral-spatial joint feature [91], [92]. Compared to the concatenation method, the authors in [94] and [95] proposed to use another SAEs for fusing the spectral and spatial features, which may further enhance the discriminative ability of spectral-spatial features.

Similar to the traditional feature extraction methods, one can also embed some prior or expected information into SAEs. Based upon the assumption that neighboring samples in the input space should have similar hidden representations, graph regularization was added to SAE for preserving such property [96], [97]. In [98], Zhou et al. imposed a local regularization via Fisher discriminant analysis on hidden layers to make the extracted features of samples from the same category close to each other while from different categories as far as possible, thus improving the discriminative ability of SAE. Meanwhile, they also added a diversity regularization term to make SAE extract compact features.

### B. Convolutional Neural Networks (CNNs)

CNNs are the most popularly adopted deep model for hyperspectral feature extraction. As shown in Fig. 9, the basic components of a CNN model include convolutional layers, pooling layers, and fully connected layers. The convolutional layers are used to extract features with convolutional kernels (filters), which can be formulated as:

$$ \mathbf{X}^l = f(\mathbf{X}^{l-1} * \mathbf{W}^l + \mathbf{b}^l), \quad (34) $$

where $\mathbf{X}^l$ is the $l$th feature maps, $\mathbf{W}^l$ and $\mathbf{b}^l$ denote the filters and biases of the $l$th layer, respectively, and '*' represents the convolutional operation. After the convolutional layer, the pooling layer is often adopted for reducing the size of the generated feature maps and producing more robust features. On the top of a CNN model, there often exist some fully connected layers, aiming at learning high-level features and outputting the final results of the network.

For hyperspectral images, CNNs can be used for extracting spectral features [99] or spatial features [100]–[102], depending on the inputs of networks. In [99], Hu et al. designed a 1-D CNN model to extract spectral features of each pixel. Compared to the traditional fully-connected networks, CNNs have weight-sharing and local-connection characteristics, making their training processes more efficient and effective. In [100], 2-D CNN was explored to extract spatial features from a local cube. Different from SAEs, CNNs do not need to reshape the cube into a vector, thus preserving as much spatial
<figure><img src="image_14.png" /><figcaption>Fig. 9. Four major categories of deep learning models.</figcaption></figure>

information as possible. However, to make full use of the representation ability of CNNs, two important issues need to be considered. The first issue is the small number of training samples but high-dimensional spectral information, which will easily lead to the overfitting problem. The second issue is the extraction of spectral-spatial joint features, which can improve the classification performance in comparison with using the spectral or spatial feature only.

For the first issue, many widely used strategies in the field of natural image classification, such as dropout and weight decay, can be naturally adopted. In addition, a lot of promising methods have been proposed in the past few years. These methods can be divided into four different classes. The first class of methods is dimensionality reduction. In [100], [101], [103], PCA was employed to extract the first principal components of hyperspectral images as inputs of CNNs, thus simplifying the network structures. Similarly, a similarity-based band selection method was used in [104]. However, these dimensionality reduction methods are independent from the following CNNs, which may lose some useful information. Different from them, Ghamisi et al. proposed a novel method to adaptively select the most informative bands suitable for the CNN model [105]. The second class of methods is data augmentation. In [101], two methods were proposed to generate virtual samples. One is to multiply a random factor and add a random noise to training samples, while the other one is to combine two given samples from the same class with proper ratios. In [106], a data augmentation method based on distance density was proposed. Recently, Kong et al. proposed a random zero setting method to generate new samples [107]. The third class of methods is transfer learning. In [108] and [109], the authors found that CNNs trained by one hyperspectral data can be transferred to another data acquired by the same sensor, and fine-tuning only a few top layers achieves satisfying results. More interestingly, the works in [110]–[112] indicated that CNNs pretrained by natural images can be directly applied to extract spatial features of hyperspectral images. The fourth class of methods is semi-supervised or even unsupervised learning. For examples, Wu and Prasad attempted to use a clustering model to obtain pseudo labels of unlabeled samples, and then combine the training samples and unlabeled samples

(with their pseudo labels) together to train their network [113].

In terms of the second issue, one popularly used method is feeding a local cube, directly cropped from the original hyperspectral image, into a CNN with 3-D convolution kernels for processing the spectral and spatial information simultaneously. The number of channels in the 3-D convolutional kernel is smaller than or equal to that of its input layer. However, the former one dramatically increases the computational complexity due to the simultaneous convolution operators in both spectral domain and spatial domain, while the latter one heavily increases the number of parameters to optimize. Another candidate method is to decouple the task of spectral-spatial feature extraction into two parts: Spectral feature extraction and spatial feature extraction. In [108] and [114], a parallel structure was employed to extract spectral-spatial features. Specifically, 1-D CNN and 2-D CNN were designed to extract spectral features and spatial features, respectively; these two features were then concatenated together and fused via a few fully-connected layers. Since 2-D CNN focuses on extracting spatial features, some redundant spectral information can be preprocessed to reduce the computational complexity. In [115], a serial structure was also used to extract spectral-spatial features. It firstly applied several 1 × 1 convolutions to extract spectral features and then fed the extracted features into several 3-D convolutions to extract spatial features.

### C. Recurrent Neural Networks (RNNs)

RNNs have been popularly employed to sequential data analysis, such as machine translation and speech recognition. Different from the feedforward neural network, RNN takes advantage of a recurrent edge for connecting the neuron to itself across time. Therefore, it is able to model the probability distribution of sequence data. To make this subsection easier to follow, we first provide a brief and general discussion on RNN. Then, we briefly describe how to use RNN specifically for the classification of hyperspectral images.

Fig. 9 shows an example of RNN. Given a sequence $x = (x_1, x_2, \cdots, x_T)$, where $x_t, t \in \{1, 2, \cdots, T\}$ generally
denotes the information at the time t, the output of the hidden layer at the tth time step is:

$$h_t = f(Ux_t + Wh_{t-1} + b_h), \quad (35)$$

where U and W represent weight matrices from the current input layer to the hidden layer and the preceding hidden layer to the current hidden layer, respectively, $h_{t-1}$ is the output of the hidden layer at the preceding time, and $b_h$ is a bias vector. According to this equation, it can be observed that the contextual relationships in the time domain are constructed via a recurrent connection. Ideally, $h_T$ will capture most of the information, and can be considered as the final feature of the sequence data. In terms of classification tasks, one often inputs $h_t$ into an output layer $o_t$, which can be described as:

$$o_t = f(Vh_t + b_o), \quad (36)$$

where $V$ is the weight matrix from the hidden layer to the output layer, and $b_o$ is a bias vector.

In recent years, RNNs have attracted more and more attention in the field of hyperspectral image feature extraction. To make full use of RNNs, one needs to first ask the following question: How to construct the sequence? An intuitive method is regarding the whole spectral bands as a sequence [116], [117]. For each pixel, its spectral values are fed into RNN from the first band to the last band, and the output of the hidden layer at the last band is the extracted spectral feature. Different from the traditional sequences in speech recognition or machine translation tasks, the succeeding bands do not depend on the preceding ones. Thus, Liu et al. also fed the spectral sequence from the last band to the first band to construct a bidirectional RNN model [116]. Another method is using a local patch or cube to construct the sequence [117]–[119]. For examples, Zhou et al. regarded the rows of each local patch, cropped from the first principal component of hyperspectral images, as a sequence, and fed them into RNN one by one to extract spatial features [117]; Zhang et al. adopted each pixel and its neighboring pixels in the cube to form a sequence [118]. These pixels were firstly sorted according to their similarities to the center pixel, and then fed into RNN sequentially to extract locally spatial features.

In real applications, the constructed sequence may be very long. Take the widely used Indian Pines data as an example, the length of the sequence is 200 (the number of spectral bands) if we use the first method mentioned above to construct the sequence. This sequence will increase the training difficulty, because of the gradients tending to either vanish or explode. In order to deal with this issue, long short-term memory (LSTM) was employed as a more sophisticated recurrent unit [116], [117], [120], [121]. The core components of LSTM are three gates: an input gate, a forget gate, and an output gate. These gates together control the flow of information in the network. Similarly, gated recurrent unit (GRU) was also employed, which has only two gates (i.e., an update gate and a reset gate). Compared to LSTM units, GRUs have less numbers of parameters, which may be more suitable for hyperspectral image feature extraction since it usually has a limited number of training samples. Another candidate scheme to address the issue is to divide the

long-term sequence into shorter sequences [121], [122]. For example, in [122], Hang et al. proposed to group the adjacent bands of hyperspectral images into subsequences, and then use RNNs to extract features from them. Since nonadjacent bands have some complementarity, they also used another RNN to fuse the extracted features.

## D. Integrated Networks

In general, AEs and RNNs are good at processing vectorized inputs, thus achieving promising results in terms of spectral feature extraction. However, both of them need to reshape the input patches or cubes into vectors during spatial feature extraction, which may destroy some spatial information. In contrast, CNNs are able to directly deal with image patches and cubes, resulting in more powerful spatial features than AEs and RNNs. It is natural to think whether we can integrate these networks together to make full use of their respective advantages. In the past few years, numerous works have been proposed in this direction.

One kind of integration method is to use each network independently, and then combining their results together [119], [121]–[123]. In [123], a parallel framework was proposed to extract spectral-spatial joint features from hyperspectral images. In this framework, SAE was employed to extract spectral features of each pixel, and CNN was used to extract spatial features from the corresponding image patch. These two results were fused by a fully connected layer. Similar to this work, Xu et al. also adopted the parallel framework but used LSTM to extract the spectral features [121]. Different from them, Hang et al. proposed a serial framework to fuse CNNs and RNNs [122]. Specifically, they used CNN to extract the spatial features from each band of hyperspectral images, and then used RNN to fuse the extracted spatial features. In [119], Shi and Pun also employed a serial framework to integrate CNN and RNN for spectral-spatial feature extraction.

Another kind of integration method is embedding the core component (i.e., convolutional operators) of CNNs into AEs or RNNs [93], [116]. In [93], an unsupervised spectral-spatial feature extraction network was proposed. The whole framework was similar to AEs, which also adopted the so-called encoder-decoder paradigm. However, the fully-connected operators in AEs were replaced by convolutional operators, so that the network can directly extract spectral-spatial joint features from cubes. In [116], Liu et al. proposed a spectral-spatial feature extraction method based on a convolutional LSTM network. Instead of fully connected operators, they also used convolutional operators in LSTM units. For a given cube, each band was fed into the convolutional LSTM unit sequentially. The convolutional operators can extract the spatial features while the recurrent operators can extract the spectral features. The whole network was optimized in an end-to-end way, thus achieving satisfactory performance.

## V. EXPERIMENTAL RESULTS

To evaluate the performance of different FE techniques, we selected four techniques from the UFE category (i.e., PCA [19], MSTV [39], OTVCA [42], LPP [48]), four techniques
from the SFE category (i.e., LDA [51], CGDA [74], LSDR [82], and JPlay [87]), and five techniques from the deep FE category (i.e., SAE [124], RNN [122], CNN [125], CAE [126], and CRNN [116]). Here, we set the tuning parameters for those algorithms before representing the experimental results.

## A. Algorithm Setup

The parameter setting usually plays a crucial role in assessing the performance of FE algorithms. Subspace dimension (or number of features, $d$) is a common parameter for all compared algorithms. Selection of the number of features is a hard task for hyperspectral image analysis. The endmember selection/extraction, subspace identification, and/or rank selection are all refer to this issue [127]–[129]. For a fair and simplified comparison, the parameter $d$ is unifiedly assigned to the same value, which is equal to the number of classes ($k$). We should note that $d$ in LDA is automatically determined as $k-1$, due to the class separability (fisher's criterion).

### 1) Unsupervised FE:

* PCA: This method is a parameter free technique.

* MSTV: In [39], all parameters are adjusted using a trial and error approach. The multi-scale parameters adjusting the degree of smoothness (as suggested in [39]) are set to 0.003, 0.02, and 0.01. The spatial scale for the structure extraction in three levels (as suggested in [39]) is set to 2, 1, and 3.

* OTVCA: This method is initialized as recommended in [42]. The tuning parameter $\lambda$ which controls the level of smoothness applied on the features is set to one percent of the maximum intensity range of the datatsets.

* LPP: The number of neighbors is set to 12. The bandwidth of the Gaussian kernel is set to 1.

### 2) Supervised FE:
A common strategy for a model selection is to run cross-validation (CV) on the training set, since the labeled samples are available in SFE. Therefore, we used the CV strategy on the following studied algorithms for parameter selection.

* LDA: This method can be viewed as a baseline for SFE. There is no additional parameter in LDA.

* CGDA: The Eq. (23) can be tuned to a regularized optimization problem, where one extra parameter – regularized $l_2$-norm – needs to be set in advance in the process of graph construction, which can be searched in the range of $\{10^{-4}, 10^{-3}, 10^{-2}, \dots, 10^1, 10^2\}$ by CV. In the experiments, 0.1 is used for all three datasets.

* LSDR: Two parameters are involved in LSDR; The standard deviation for the Gaussian function and the regularization parameter which are selected in the range of $\{0.05, 0.1, \dots, 0.95, 1\}$ and $\{10^{-2}, 10^{-1}, \dots, 10^1, 10^2\}$, respectively, using CV. Finally, $\sigma$ and $\lambda$ are both set to 1 in our experiments.

* JPlay: There are three regularization parameters ($\alpha, \beta$, and $\gamma$) that need to be set in the JPlay model (32). With the CV conducted on the training set of three different datasets, the regularization parameters are selected in the ranges of $\{10^{-2}, 10^{-1}, \dots, 10^1, 10^2\}$, yielding the final setting of $(\alpha, \beta, \gamma) = (0.1, 1, 1)$ for the first dataset,

($\alpha, \beta, \gamma$) = (0.1, 0.1, 1) for the second dataset, and
($\alpha, \beta, \gamma$) = (0.1, 1, 1) for the last dataset.

### 3) Deep FE:

* SAE: The input of SAE is the original spectral information of each pixel. Three hidden layers are used. The numbers of neurons from the first to the third hidden layer are set to 32, 64, and 128, respectively. ReLU is adopted as the activation function for each hidden layer.

* CNN: The input of CNN is a small cube with a size of $16 \times 16 \times p$, where $p$ represents the number of spectral bands for each data. Three convolutional layers are used. Each convolutional layer is sequentially followed by a batch normalization layer, a ReLU activation function, and a max-pooling layer. Note that the last pooling layer is an adaptive max-pooling layer, making the output size equals to $1 \times 1$ for any input sizes. The kernel size for each convolution is $3 \times 3$, and the numbers of kernels from the first to the third convolutions are set to 32, 64, and 128, respectively. Padding operators are used to preserve the spatial size after each convolutional operator.

* PCNN: PCA is applied prior to CNN to reduce the spectral dimension of the HSI. The number of reduced dimension by PCA is set to the number of classes ($k$). The input cube for CNN is of size $16 \times 16 \times k$.

* RNN: The input of RNN is the same as the input of SAE. Two recurrent layers with GRU are employed. The number of neurons in each recurrent layer is set to 128.

* Integrated Networks: Convolutional AE (CAE) and convolutional RNN (CRNN) are selected as two representative integrated networks. The input for them is the same as that for the CNN. For CAE, three convolutional layers and three de-convolutional layers are adopted. All of them use $3 \times 3$ kernels. The numbers of kernels from the first to the third convolutional layers are set to 32, 64, and 128, respectively. In contrast, the numbers of kernels from the first and the third de-convolutional layers are set to 64, 32, and $p$, respectively. Similar to [116], CRNN adopts two recurrent layers with convolutional LSTM units. For both recurrent layers, $3 \times 3$ convolutional kernels are applied. The numbers of kernels for the first and the second recurrent layers are set to 32 and 64, respectively.

All of the above deep learning related models are implemented in the PyTorch framework. In order to optimize them, we use the Adam algorithm with default parameters. The batch size, the learning rate, and the number of training epochs are set to 128, 0.001, and 200, respectively. To reduce the effects of random initialization, all of the deep learning models are repeated 5 times, and the mean values are reported.

### 4) Random Forest Classifier:
Apart from the deep FE techniques, all the other FE techniques use random forest (RF) to perform the classification task. The number of trees selected for RF is set to 200. We approximately set the number of the prediction variable to the square root of the number of input bands.

## B. Performance of the FE techniques on three HSIs

We have applied the FE techniques on the three hyperspectral datasets, i.e., Indian Pinse 2010, Houston 2013, and
<figure><img src="image_15.png" /><figcaption>Fig. 10. Classification maps obtained on the extracted features from Indian Pines 2010 data. From each category the method with highest OA is shown for the demonstration along with the one obtained from the spectral bands.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 11. Classification maps obtained on the extracted features from Houston University 2013 data. From each category the method with highest OA is shown for the demonstration along with the one obtained from the spectral bands.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 12. Classification maps obtained on the extracted features from Houston University 2018 data. From each category the method with highest OA is shown for the demonstration along with the one obtained from the spectral bands.</figcaption></figure>
<table><caption>TABLE V
CLASSIFICATION ACCURACIES OBTAINED ON FEATURES EXTRACTED FROM THE INDIAN PINES 2010 DATASET USING DIFFERENT SHALLOW AND DEEP FEATURE EXTRACTION TECHNIQUES.</caption><thead><tr><th colspan="15">Indian Pines 2010</th></tr><tr><th rowspan="3"></th><th rowspan="3">Spectral</th><th colspan="8">Shallow FE</th><th colspan="6">Deep FE</th></tr><tr><th colspan="4">UFE</th><th colspan="4">SFE</th><th colspan="6">SFE</th></tr><tr><th>PCA</th><th>MSTV</th><th>OTVCA</th><th>LPP</th><th>LDA</th><th>CGDA</th><th>LSDR</th><th>JPlay</th><th>SAE</th><th>RNN</th><th>CNN</th><th>CAE</th><th>CRNN</th><th>PCNN</th></tr></thead><tbody><tr><td>1</td><td>0.9260</td><td>0.9064</td><td><strong>0.9992</strong></td><td>0.9260</td><td>0.8850</td><td>0.9628</td><td>0.8144</td><td>0.8459</td><td>0.9230</td><td>0.9327</td><td>0.8829</td><td>0.9275</td><td>0.9432</td><td>0.9590</td><td>0.9397</td></tr><tr><td>2</td><td>0.8769</td><td>0.9976</td><td><strong>1.0000</strong></td><td>0.9976</td><td><strong>1.0000</strong></td><td>0.9961</td><td>0.8933</td><td>0.9984</td><td>0.9573</td><td>0.9178</td><td>0.9544</td><td>0.9961</td><td>0.8596</td><td>0.9998</td></tr><tr><td>3</td><td>0.8862</td><td>0.9724</td><td>0.9724</td><td>0.9897</td><td>0.9724</td><td>0.9724</td><td>0.9759</td><td>0.9655</td><td>0.9683</td><td>0.9405</td><td>0.9890</td><td><strong>0.9986</strong></td><td>0.8241</td><td>0.9938</td></tr><tr><td>4</td><td>0.6888</td><td>0.7762</td><td><strong>1.0000</strong></td><td>0.8953</td><td>0.8742</td><td>0.9270</td><td>0.7474</td><td>0.8194</td><td>0.8300</td><td>0.7508</td><td>0.7621</td><td>0.8840</td><td>0.8822</td><td>0.8569</td><td>0.8930</td></tr><tr><td>5</td><td>0.8058</td><td><strong>0.8855</strong></td><td>0.8039</td><td>0.8151</td><td>0.8682</td><td>0.8802</td><td>0.8096</td><td>0.8394</td><td>0.8665</td><td>0.8488</td><td>0.8474</td><td>0.8692</td><td>0.8570</td><td>0.8638</td><td>0.8706</td></tr><tr><td>6</td><td>0.8172</td><td>0.8797</td><td><strong>0.9946</strong></td><td>0.7094</td><td>0.9739</td><td>0.7883</td><td>0.8284</td><td>0.9397</td><td>0.9418</td><td>0.9013</td><td>0.9204</td><td>0.9268</td><td>0.9167</td><td>0.6975</td><td>0.9127</td></tr><tr><td>7</td><td>0.4170</td><td>0.5954</td><td>0.6792</td><td><strong>0.7166</strong></td><td>0.6985</td><td><strong>0.7166</strong></td><td>0.6845</td><td>0.7030</td><td>0.6958</td><td>0.6265</td><td>0.5795</td><td>0.6507</td><td>0.6818</td><td>0.6348</td><td>0.7103</td></tr><tr><td>8</td><td>0.2530</td><td>0.2583</td><td>0.2599</td><td>0.2768</td><td>0.2961</td><td>0.2955</td><td>0.2431</td><td>0.2952</td><td>0.2758</td><td>0.5349</td><td>0.2840</td><td>0.8776</td><td>0.6776</td><td><strong>0.9934</strong></td><td>0.4725</td></tr><tr><td>9</td><td>0.6545</td><td>0.7498</td><td>0.7048</td><td>0.7943</td><td><strong>0.8913</strong></td><td>0.8457</td><td>0.7971</td><td>0.8419</td><td>0.8142</td><td>0.8732</td><td>0.8533</td><td>0.8302</td><td>0.8336</td><td>0.8194</td><td>0.8621</td></tr><tr><td>10</td><td>0.8229</td><td>0.9406</td><td>0.9594</td><td>0.9368</td><td>0.9019</td><td><strong>0.9804</strong></td><td>0.8514</td><td>0.7761</td><td>0.9663</td><td>0.9015</td><td>0.8096</td><td>0.8368</td><td>0.8752</td><td>0.7946</td><td>0.9289</td></tr><tr><td>11</td><td>0.6658</td><td>0.8402</td><td>0.9224</td><td><strong>0.9945</strong></td><td>0.8943</td><td>0.9288</td><td>0.7195</td><td>0.7651</td><td>0.8052</td><td>0.8121</td><td>0.7414</td><td>0.7440</td><td>0.7633</td><td>0.6492</td><td>0.9165</td></tr><tr><td>12</td><td>0.9985</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9995</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9998</td><td>0.9765</td><td>0.9945</td><td>0.9838</td><td>0.9748</td><td>0.9991</td></tr><tr><td>13</td><td>0.9468</td><td><strong>0.9962</strong></td><td>0.9879</td><td>0.9888</td><td>0.9925</td><td>0.9830</td><td>0.9580</td><td>0.9738</td><td>0.9819</td><td>0.9621</td><td>0.9427</td><td>0.9925</td><td>0.9930</td><td>0.9892</td><td>0.9959</td></tr><tr><td>14</td><td>0.8783</td><td>0.9000</td><td>0.9615</td><td>0.9145</td><td>0.9344</td><td>0.9174</td><td>0.8756</td><td>0.8628</td><td>0.8953</td><td>0.9094</td><td>0.9030</td><td>0.9984</td><td>0.9985</td><td>0.8981</td><td><strong>0.9993</strong></td></tr><tr><td>15</td><td>0.9333</td><td>0.9667</td><td>0.9511</td><td>0.9933</td><td>0.9489</td><td>0.9756</td><td>0.9333</td><td>0.9311</td><td>0.9600</td><td>0.9307</td><td>0.8119</td><td>0.9947</td><td>0.9942</td><td>0.7556</td><td><strong>0.9978</strong></td></tr><tr><td>16</td><td>0.3735</td><td>0.2036</td><td>0.2885</td><td>0.1601</td><td>0.5217</td><td>0.1719</td><td>0.3439</td><td>0.4901</td><td>0.4466</td><td>0.2053</td><td>0.2060</td><td>0.0980</td><td>0.1700</td><td><strong>0.6028</strong></td><td>0.1051</td></tr><tr><td>AA</td><td>0.7465</td><td>0.8043</td><td>0.8428</td><td>0.8194</td><td><strong>0.8532</strong></td><td>0.8341</td><td>0.7861</td><td>0.8089</td><td>0.8367</td><td>0.8197</td><td>0.7737</td><td>0.8480</td><td>0.8478</td><td>0.8233</td><td>0.8498</td></tr><tr><td>OA</td><td>0.7868</td><td>0.8598</td><td>0.8561</td><td>0.8378</td><td><strong>0.9112</strong></td><td>0.8748</td><td>0.8370</td><td>0.8748</td><td>0.8829</td><td>0.8836</td><td>0.8655</td><td>0.8945</td><td>0.8911</td><td>0.8525</td><td>0.9018</td></tr><tr><td>κ</td><td>0.7390</td><td>0.8297</td><td>0.8247</td><td>0.8054</td><td><strong>0.8909</strong></td><td>0.8481</td><td>0.8010</td><td>0.8466</td><td>0.8571</td><td>0.8580</td><td>0.8355</td><td>0.8716</td><td>0.8673</td><td>0.8213</td><td>0.8802</td></tr></tbody></table>

<table><caption>TABLE VI
CLASSIFICATION ACCURACIES OBTAINED ON FEATURES EXTRACTED FROM THE HOUSTON UNIVERSITY 2013 DATASET USING DIFFERENT SHALLOW AND DEEP FEATURE EXTRACTION TECHNIQUES.</caption><thead><tr><th colspan="15">Houston 2013</th></tr><tr><th rowspan="3"></th><th rowspan="3">Spectral</th><th colspan="8">Shallow FE</th><th colspan="6">Deep FE</th></tr><tr><th colspan="4">UFE</th><th colspan="4">SFE</th><th colspan="6">SFE</th></tr><tr><th>PCA</th><th>MSTV</th><th>OTVCA</th><th>LPP</th><th>LDA</th><th>CGDA</th><th>LSDR</th><th>JPlay</th><th>SAE</th><th>RNN</th><th>CNN</th><th>CAE</th><th>CRNN</th><th>PCNN</th></tr></thead><tbody><tr><td>1</td><td>0.8262</td><td><strong>0.8272</strong></td><td>0.8025</td><td>0.8205</td><td>0.8110</td><td>0.8177</td><td>0.8139</td><td>0.8120</td><td>0.7768</td><td>0.8217</td><td>0.8182</td><td>0.8104</td><td>0.8154</td><td>0.8245</td><td>0.8089</td></tr><tr><td>2</td><td>0.8318</td><td>0.8393</td><td>0.8412</td><td>0.8515</td><td>0.8214</td><td>0.8355</td><td>0.8327</td><td>0.8553</td><td><strong>0.9662</strong></td><td>0.8274</td><td>0.8153</td><td>0.8425</td><td>0.8167</td><td>0.8412</td><td>0.8293</td></tr><tr><td>3</td><td>0.9782</td><td><strong>1.0000</strong></td><td>0.9822</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9980</td><td>0.9895</td><td>0.9939</td><td>0.8594</td><td>0.7731</td><td>0.9156</td><td>0.8432</td></tr><tr><td>4</td><td>0.9138</td><td>0.9081</td><td>0.7633</td><td>0.8873</td><td>0.9479</td><td>0.8920</td><td>0.9053</td><td>0.8864</td><td>0.9564</td><td><strong>0.9773</strong></td><td>0.9040</td><td>0.9170</td><td>0.9153</td><td>0.9129</td><td>0.9159</td></tr><tr><td>5</td><td>0.9659</td><td>0.9886</td><td>0.9915</td><td><strong>0.9991</strong></td><td>0.9867</td><td>0.9384</td><td>0.9915</td><td>0.9688</td><td>0.9782</td><td>0.9438</td><td>0.9389</td><td>0.9699</td><td>0.9585</td><td>0.9881</td><td>0.9824</td></tr><tr><td>6</td><td>0.9930</td><td>0.9930</td><td>0.9580</td><td>0.9580</td><td>0.9790</td><td><strong>1.0000</strong></td><td>0.8741</td><td>0.9860</td><td>0.9930</td><td>0.9874</td><td>0.9678</td><td>0.8769</td><td>0.9776</td><td>0.9483</td><td>0.9497</td></tr><tr><td>7</td><td>0.7463</td><td>0.8927</td><td>0.6362</td><td>0.7090</td><td><strong>0.9123</strong></td><td>0.7901</td><td>0.8535</td><td>0.8526</td><td>0.7817</td><td>0.7293</td><td>0.7392</td><td>0.8802</td><td>0.8694</td><td>0.8642</td><td>0.8627</td></tr><tr><td>8</td><td>0.3305</td><td>0.4606</td><td>0.5992</td><td>0.6724</td><td>0.4311</td><td>0.7379</td><td>0.4302</td><td>0.4710</td><td>0.7806</td><td>0.3792</td><td>0.4153</td><td>0.6344</td><td>0.6762</td><td>0.5305</td><td><strong>0.8351</strong></td></tr><tr><td>9</td><td>0.6771</td><td>0.7885</td><td>0.8706</td><td><strong>0.9008</strong></td><td>0.7413</td><td>0.6449</td><td>0.7186</td><td>0.6752</td><td>0.7592</td><td>0.7145</td><td>0.7367</td><td>0.8595</td><td>0.8540</td><td>0.8404</td><td>0.8691</td></tr><tr><td>10</td><td>0.4295</td><td>0.4749</td><td>0.6612</td><td><strong>0.8398</strong></td><td>0.4595</td><td>0.4662</td><td>0.4826</td><td>0.6014</td><td>0.5556</td><td>0.5373</td><td>0.5674</td><td>0.5782</td><td>0.4514</td><td>0.6168</td><td>0.7913</td></tr><tr><td>11</td><td>0.7011</td><td>0.7268</td><td>0.9820</td><td><strong>0.9924</strong></td><td>0.7306</td><td>0.7239</td><td>0.7287</td><td>0.5806</td><td>0.6983</td><td>0.6231</td><td>0.7250</td><td>0.7417</td><td>0.7292</td><td>0.6186</td><td>0.7913</td></tr><tr><td>12</td><td>0.5485</td><td>0.9145</td><td>0.7349</td><td><strong>0.9625</strong></td><td>0.7560</td><td>0.6513</td><td>0.7656</td><td>0.5687</td><td>0.7858</td><td>0.6305</td><td>0.7606</td><td>0.9379</td><td>0.9402</td><td>0.8440</td><td>0.9593</td></tr><tr><td>13</td><td>0.6140</td><td>0.7754</td><td>0.6982</td><td>0.7789</td><td>0.8105</td><td>0.6105</td><td>0.7719</td><td>0.6702</td><td>0.7509</td><td>0.4516</td><td>0.6656</td><td>0.8835</td><td><strong>0.8968</strong></td><td>0.8414</td><td>0.8765</td></tr><tr><td>14</td><td>0.9838</td><td>0.9919</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9960</td><td>0.9919</td><td>0.9879</td><td>0.9595</td><td>0.9879</td><td>0.9692</td><td>0.9850</td><td>0.9943</td><td>0.9773</td><td>0.9603</td><td>0.9968</td></tr><tr><td>15</td><td>0.9789</td><td>0.9746</td><td><strong>1.0000</strong></td><td>0.9789</td><td>0.9746</td><td>0.9831</td><td>0.9852</td><td>0.9514</td><td>0.9831</td><td>0.9732</td><td>0.9607</td><td>0.8072</td><td>0.7471</td><td>0.9345</td><td>0.8592</td></tr><tr><td>AA</td><td>0.7679</td><td>0.8371</td><td>0.8347</td><td><strong>0.8901</strong></td><td>0.8239</td><td>0.8056</td><td>0.8094</td><td>0.7878</td><td>0.8532</td><td>0.7716</td><td>0.7976</td><td>0.8388</td><td>0.8350</td><td>0.8210</td><td>0.8664</td></tr><tr><td>OA</td><td>0.7278</td><td>0.8058</td><td>0.8088</td><td><strong>0.8753</strong></td><td>0.7874</td><td>0.7745</td><td>0.7789</td><td>0.7524</td><td>0.8280</td><td>0.7436</td><td>0.7646</td><td>0.8239</td><td>0.8184</td><td>0.7921</td><td>0.8526</td></tr><tr><td>κ</td><td>0.7076</td><td>0.7895</td><td>0.7923</td><td><strong>0.8648</strong></td><td>0.7700</td><td>0.7552</td><td>0.7604</td><td>0.7315</td><td>0.8134</td><td>0.7235</td><td>0.7469</td><td>0.8096</td><td>0.8036</td><td>0.7761</td><td>0.8404</td></tr></tbody></table>

Houston 2018, and the classification accuracies including class accuracies, AA, OA, and $\kappa$ coefficient are all shown in Tables V, VI, and VII, respectively. The results are first discussed within the categories and then between different categories. We should note that the results and discussions are in terms of classification accuracies obtained from the classification of the HSIs.

### 1) Unsupervised FE:

* PCA demonstrates the poorest performance compared with the other techniques, however, it considerably improves the classification accuracies compared with the

results obtained by applying the RF on the spectral bands. One of the main disadvantages of the PCA is that, it does not take into account the noise and, therefore, the extracted features with having lower variance are often degraded by different types of noise existing in the hyperspectral image [130]. Additionally, PCA only takes into account the spectral correlation and it entirely neglects the spatial (neighboring) information.

* LPP considerably outperforms the other UFE techniques for the Indian Pines dataset. However, in the case of Houston datasets, it provides very poor results. LPP
<table><caption>TABLE VII
CLASSIFICATION ACCURACIES OBTAINED ON FEATURES EXTRACTED FROM THE HOUSTON UNIVERSITY 2018 DATASET USING DIFFERENT SHALLOW AND DEEP FEATURE EXTRACTION TECHNIQUES.</caption><thead><tr><th rowspan="3"></th><th colspan="12">Houston 2018</th></tr><tr><th colspan="12">Shallow FE</th><th colspan="6">Deep FE</th></tr><tr><th colspan="3">UFE</th><th colspan="3">SFE</th><th colspan="6">SFE</th></tr><tr><th></th><th>Spectral</th><th>PCA</th><th>MSTV</th><th>OTVCA</th><th>LPP</th><th>LDA</th><th>CGDA</th><th>LSDR</th><th>JPlay</th><th>SAE</th><th>RNN</th><th>CNN</th><th>CAE</th><th>CRNN</th><th>PCNN</th></tr></thead><tbody><tr><td>1</td><td>0.3088</td><td><strong>0.8781</strong></td><td>0.0536</td><td>0.6842</td><td>0.6618</td><td>0.6256</td><td>0.7575</td><td>0.7969</td><td>0.5991</td><td>0.7940</td><td>0.5702</td><td>0.7516</td><td>0.4428</td><td>0.6338</td><td>0.6638</td></tr><tr><td>2</td><td>0.7603</td><td>0.8396</td><td>0.7046</td><td>0.6376</td><td>0.8122</td><td>0.8474</td><td>0.8076</td><td>0.7747</td><td>0.8347</td><td>0.7893</td><td>0.6975</td><td>0.8173</td><td><strong>0.8849</strong></td><td>0.8707</td><td>0.8376</td></tr><tr><td>3</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9972</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9972</td><td>0.7739</td><td>0.8482</td><td>0.9924</td><td>0.8045</td></tr><tr><td>4</td><td>0.9134</td><td>0.9494</td><td>0.6238</td><td>0.6775</td><td>0.9453</td><td>0.9059</td><td>0.9265</td><td>0.9276</td><td>0.9298</td><td>0.9221</td><td>0.8613</td><td>0.9444</td><td>0.9362</td><td>0.9439</td><td><strong>0.9595</strong></td></tr><tr><td>5</td><td>0.4119</td><td>0.4668</td><td>0.2676</td><td>0.1679</td><td>0.4728</td><td>0.5258</td><td>0.4661</td><td>0.4289</td><td>0.3971</td><td>0.4982</td><td>0.4040</td><td>0.4330</td><td>0.5396</td><td><strong>0.5404</strong></td><td>0.4800</td></tr><tr><td>6</td><td>0.2570</td><td>0.2990</td><td><strong>0.3835</strong></td><td>0.3164</td><td>0.3008</td><td>0.2910</td><td>0.2776</td><td>0.2726</td><td>0.2780</td><td>0.2585</td><td>0.2537</td><td>0.3050</td><td>0.3080</td><td>0.2902</td><td>0.3377</td></tr><tr><td>7</td><td>0.3025</td><td>0.3025</td><td>0.3025</td><td>0.3025</td><td>0.3025</td><td>0.3025</td><td>0.3109</td><td>0.3025</td><td>0.2857</td><td>0.3025</td><td>0.3025</td><td>0.2908</td><td>0.2723</td><td>0.2997</td><td><strong>0.3176</strong></td></tr><tr><td>8</td><td>0.7657</td><td>0.7675</td><td>0.7599</td><td>0.7417</td><td>0.7785</td><td>0.7849</td><td>0.7544</td><td>0.7518</td><td>0.7771</td><td>0.7216</td><td>0.7356</td><td>0.8538</td><td>0.8583</td><td>0.8092</td><td><strong>0.8677</strong></td></tr><tr><td>9</td><td>0.3849</td><td>0.3877</td><td>0.5767</td><td>0.5990</td><td>0.4887</td><td>0.3917</td><td>0.3672</td><td>0.5255</td><td>0.5877</td><td>0.6302</td><td>0.4186</td><td>0.7970</td><td>0.7371</td><td>0.3717</td><td><strong>0.8659</strong></td></tr><tr><td>10</td><td>0.3603</td><td>0.4360</td><td>0.3747</td><td>0.4491</td><td>0.4230</td><td>0.4086</td><td>0.3790</td><td>0.3497</td><td>0.4010</td><td>0.3819</td><td>0.3465</td><td><strong>0.5902</strong></td><td>0.4957</td><td>0.5484</td><td>0.5778</td></tr><tr><td>11</td><td>0.4162</td><td>0.4792</td><td><strong>0.7862</strong></td><td>0.7596</td><td>0.5085</td><td>0.4667</td><td>0.4266</td><td>0.4422</td><td>0.5359</td><td>0.4143</td><td>0.4699</td><td>0.5456</td><td>0.5781</td><td>0.6048</td><td>0.5948</td></tr><tr><td>12</td><td>0.0132</td><td>0.0046</td><td>0.0093</td><td>0.0077</td><td>0.0070</td><td>0.0023</td><td>0.0170</td><td>0.0000</td><td>0.0302</td><td>0.0152</td><td>0.0697</td><td>0.0511</td><td>0.0477</td><td><strong>0.0927</strong></td><td>0.0579</td></tr><tr><td>13</td><td>0.4525</td><td>0.5556</td><td>0.4238</td><td>0.4090</td><td>0.5442</td><td>0.5164</td><td>0.5707</td><td>0.5603</td><td>0.5324</td><td>0.4523</td><td>0.4789</td><td>0.5148</td><td>0.5619</td><td>0.4246</td><td><strong>0.5811</strong></td></tr><tr><td>14</td><td>0.3019</td><td>0.2629</td><td>0.5460</td><td>0.4665</td><td>0.3651</td><td>0.4152</td><td>0.2294</td><td>0.2073</td><td>0.3212</td><td>0.3789</td><td>0.3309</td><td>0.5289</td><td><strong>0.6763</strong></td><td>0.3375</td><td>0.5705</td></tr><tr><td>15</td><td>0.6303</td><td>0.4721</td><td>0.4457</td><td>0.4887</td><td>0.4602</td><td>0.5549</td><td>0.4180</td><td>0.5234</td><td>0.5944</td><td>0.5197</td><td>0.5289</td><td>0.6277</td><td>0.6476</td><td>0.6447</td><td><strong>0.6591</strong></td></tr><tr><td>16</td><td>0.6412</td><td>0.7611</td><td>0.6220</td><td>0.7648</td><td>0.7559</td><td>0.5888</td><td>0.7374</td><td>0.6403</td><td>0.6688</td><td>0.7457</td><td>0.7086</td><td>0.8498</td><td>0.7594</td><td>0.7173</td><td><strong>0.8572</strong></td></tr><tr><td>17</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9545</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td><td>0.9545</td><td>0.8909</td><td><strong>1.0000</strong></td><td><strong>1.0000</strong></td></tr><tr><td>18</td><td>0.4983</td><td>0.6885</td><td>0.6576</td><td>0.5197</td><td>0.7140</td><td>0.6581</td><td>0.6625</td><td>0.5686</td><td>0.7366</td><td>0.5346</td><td>0.6080</td><td>0.6365</td><td>0.5981</td><td><strong>0.7692</strong></td><td>0.7020</td></tr><tr><td>19</td><td>0.5265</td><td>0.6363</td><td>0.9060</td><td>0.8777</td><td>0.7323</td><td>0.6989</td><td>0.6100</td><td>0.6266</td><td>0.6771</td><td>0.6569</td><td>0.6545</td><td>0.9102</td><td>0.8960</td><td>0.8006</td><td><strong>0.9476</strong></td></tr><tr><td>20</td><td>0.4444</td><td>0.8904</td><td><strong>0.9706</strong></td><td>0.5253</td><td>0.8479</td><td>0.9189</td><td>0.6797</td><td>0.5955</td><td>0.5393</td><td>0.5388</td><td>0.4801</td><td>0.6246</td><td>0.5566</td><td>0.4879</td><td>0.6519</td></tr><tr><td>AA</td><td>0.5195</td><td>0.6039</td><td>0.5707</td><td>0.5696</td><td>0.6060</td><td>0.5952</td><td>0.5676</td><td>0.5647</td><td>0.5863</td><td>0.5777</td><td>0.5458</td><td>0.6400</td><td>0.6268</td><td>0.6090</td><td><strong>0.6667</strong></td></tr><tr><td>OA</td><td>0.4634</td><td>0.5101</td><td>0.5750</td><td>0.5899</td><td>0.5552</td><td>0.5027</td><td>0.4825</td><td>0.5492</td><td>0.5944</td><td>0.5938</td><td>0.4851</td><td>0.7278</td><td>0.6969</td><td>0.5116</td><td><strong>0.7728</strong></td></tr><tr><td>κ</td><td>0.3732</td><td>0.4317</td><td>0.4833</td><td>0.4974</td><td>0.4714</td><td>0.4231</td><td>0.4018</td><td>0.4560</td><td>0.5037</td><td>0.4948</td><td>0.3936</td><td>0.6474</td><td>0.6124</td><td>0.4372</td><td><strong>0.7011</strong></td></tr></tbody></table>

incorporates the spatial information using the manifold learning process and by constructing the neighboring graph [48].

▷ OTVCA outperforms the other UFE techniques for Houston datasets. In the case of Houston 2013, the improvements are very considerable. OTVCA is robust to noise due to the signal model which takes into account the noise and model's errors. Additionally, OTVCA exploits the spatial correlation by incorporating the TV penalty and, therefore, the extracted features are piece-wise smooth and have high SNR [42]. Overall, it can be observed that OTVCA, which is a candidate from the low-rank reconstruction techniques, generally provides better classification accuracies than other the UFE techniques.

## 2) Supervised FE:

▷ LDA versus Spectral Classifier (RF): With the embedding of supervised information, LDA obviously performs better than the situation where RF is directly applied to the spectral signatures in terms of the overall performance and individual accuracies for most materials. This indicates the effectiveness of SFE to a great extent.

▷ LDA versus CGDA: Although the classification performance of CGDA is inferior to that of LDA from the whole perspective, yet the advantage of CGDA mainly lies in its automation in computing the similarity (or connectivity) between samples. This could lead to a relatively stable FE, particularly in large-scale and more complex hyperspectral scenes. Due to the data-driven graph embedding, CGDA yields a lower running speed than LDA in the process of model training.

▷ LDA versus LSDR: Intuitively, LSDR provides competitive classification performance with LDA. However, LSDR is time consuming due to the distribution matching between input samples and labels. The requirement to estimate the statistical distribution also limits the LSDR's stability, especially when the training set is available on a small scale (e.g., for the Indian Pines 2010 and Houston 2013 datasets).

▷ LDA versus JPlay: Unlike the conventional regression techniques, JPlay is capable of extracting semantically meaningful and robust features, due to the multi-layered structure and the self-reconstruction constraint (32). Quantitatively speaking, JPlay outperforms the other SFE methods. The CV provides a feasible solution to automatically determine the parameter combination in JPlay. Despite the ADMM solver designed for speeding up the optimization process such a multi-layered parameter update inevitably suffers from high computational cost.

## 3) Deep FE:

▷ Spectral versus Spectral-Spatial models: Most of the spectral-spatial models (i.e., CNN, PCNN, CAE, and CRNN) achieve superior performance than spectral models (i.e., SAE and RNN) in terms of AA, OA, and Kappa due to the joint use of spectral and spatial information. This indicates that, besides the rich spectral information, spatial information is also important for hyperspectral image classification.

▷ PCNN and CNN versus CAE and CRNN: Similar to SAE, CAE focuses on image reconstruction rather than classification. In contrast, PCNN and CNN are exclu-
sively designed for the classification task, so they are able to learn more discriminative features than CAE, leading to better classification performance especially on the Houston 2018 dataset. Although CRNN also focuses on the classification task, it has a higher number of parameters to train. Using the same number of training samples and epochs, PCNN and CNN can achieve better results than CRNN in terms of AA, OA, and Kappa.

▷ **PCNN versus CNN:** PCNN outperforms CNN in terms of the classification accuracies for all three datasets. We should note that the improvements are substantial in the case of the Houston 2013 and 2018 datasets. Due to the use of PCA, most of the redundant spectral information is reduced. Therefore, the number of trainable parameters in PCNN is smaller than that of CNN, making it easier to learn under the same condition.

4) **Shallow UFE Versus Shallow SFE:** For all three datasets used in the experiments, the UFE techniques provide better classification accuracies than the SFE techniques. Unlike the SFE, the UFE tends to pay more attentions on spatial-spectral information extraction, owing to fully considering all samples of HSI as the model input. Conversely, the performance of SFE is, to a great extent, limited by the ability to largely gather HSI ground sampling. A direct evident is given in Tables V, VI, and VII. For the Indian Pines 2010 and Houston 2018 datasets where more training samples are available, SFE-based methods hold the competitive results with UFE-based ones, while for the Houston 2013 dataset, the classification performance of SFE is relatively inferior to that of UFE, due to the small-scale training set. Considering the low number of ground samples often available in HSI applications the experimental results confirm the advantage of UFE over SFE for HSI feature extraction.

5) **Shallow FE versus Deep FE:** At the first glance, the shallow FE approaches slightly outperform the deep FE techniques for the two datasets, i.e., Indian Pines 2010 and Houston 2013. However, a deep comparison reveals that some deep FE techniques, such as CNN-based FE, provide consistency and good performance over all three datasets. Additionally, when the dimension-reduced step (e.g., using PCA) is applied prior to the CNN technique, the resulting PCNN yields by far the second highest accuracies in the case of Indian Pines 2010 and Houston 2013 datasets (only moderately lower than LPP or OTVCA, respectively), and simultaneously obtains the best performance on the Houston 2018 dataset. It is worth mentioning that CNN-based FE methods obtain at least 10% increase over the shallow ones in the case of Houston 2018. This could be due to the high nonlinear behaviour of this dataset which contains 20 land cover classes. The main factors for CNN-based FE methods to obtain around 20% of improvement over those shallow FE methods on the Houston 2018 dataset are the availability of the sufficient training samples and modeling the spatial information of HSI well.

6) **Comparisons of the Land Cover Maps:** Figs. 10, 11, and 12 compare the classification maps for Indian Pines, Houston University 2013, and Houston University 2018, respectively. The figures compare the maps obtained from methods which provide the highest OA from each category (i.e., Shallow UFE,

shallow SFE, and deep FE) along with the map obtained from the spectral classifier (HSI). Additionally, we depicted the maps obtained by CNN for all three datasets since it provides the highest OA among the deep FE techniques which do not exploit a reduction step.

Overall, the classification maps of either unsupervised or supervised FE-based approaches (e.g., LPP, JPlay, CNN, PCNN) are smoother compared to the results only using HSI that tends to generate sparse mislabeled pixels. More specifically, the classification maps generated by spectral-spatial FE-based methods, e.g., OTVCA, CNN, and PCNN, are usually a bit over-smoothed, leading to the creation of fake structures, especially for the Indian Pines 2010 and Houston 2018 datasets. In the case of OTVCA, the over smoothing can be avoided by decreasing the tuning parameter. In contrast, JPlay obtains relatively desirable classification maps, despite the lack of spatial information modeling. It is worth mentioning that the JPlay algorithm can maintain the structural information for Houston 2013 in the shadow covered region where pixels at some bands are considerably attenuated. This is due to the elimination of the spectral variability using the self-reconstruction regularization (the third term in (32)) and the multi-layered linearized regression technique.

### C. Performance with respect to the Number of Training Samples

In this subsection, we investigate the performance of the FE techniques in terms of classification accuracies with respect to the number of training samples. As we have already stated, this analysis is of great interest due to two main reasons; first, ground sample acquisition and measurements are often cumbersome and could be impossible in cases for which the target area is not reachable. Additionally, the limited number of samples not only affects the performance of the supervised classifiers but also the supervised features extraction techniques since they are highly reliable on the number of training samples. Therefore, in this experiment, we perform an analysis on the Houston University dataset 2018 by comparing the performances of the feature extraction techniques by selecting 10, 25, 50, and 100 training samples randomly. Fig 13 compares the OAs obtained by applying RF on the spectral bands (labeled by HSI), the features extracted by OTVCA, and JPlay along with the OAs obtained by CNN, and PCNN. The results are mean values over 10 experiments by selecting the samples randomly (the standard deviations are shown by the error bars). The outcomes of the experiment can be summarized as follows:

▷ The SFE technique (i.e., JPlay) improves the OAs compared to the spectral classifier. However, it provides much lower OA compared with the UFE and deep FE for all cases. Two aspects might explain this point. One is that JPlay fails to model spatial and contextual information, and another is that although JPlay attempts to enhance the reorientation ability of the features via multi-layered linear mapping, yet it is still incomparable to the nonlinear deep FE-based techniques, particularly when the number of samples are increased.
<figure><img src="image_18.png" /><figcaption>Fig. 13. The classification accuracies w.r.t. the number of the training samples on Houston University 2018 data. The results shown are means over 10 experiments and the standard deviations are shown by the error bars.</figcaption></figure>

▷ In this experiment, the UFE technique (i.e. OTVCA) and the deep FE one, CNN, are performed similarly in terms of classification accuracies. Compared with the results given in Table VII, it can be observed that the random selection of the training samples over the entire class regions from the ground reference considerably improves the performance of RF applied on the features extracted by OTVCA. This is often due to the lack of a parameter selection technique to select the optimum parameter for the OTVCA algorithm which could lead to over-smoothing on the features.

▷ The deep learning technique (i.e., PCNN) after using the reduction (i.e., PCA) provides very high OA for all the cases. Comparing the results with CNN (i.e., without using the PCA reduction) confirms the advantage of using the reduction stage prior to deep learning techniques.

# VI. CONCLUSION AND SUMMERY

In the past decade, hyperspectral image FE has considerably evolved leading to three main research lines (i.e., shallow unsupervised, shallow supervised, and deep feature extraction approaches) that include the majority of feature extraction techniques, which were studied in this paper. The paper has systematically provided a technical overview of the state-of-the-art techniques proposed in the literature by categorizing the aforementioned three focuses into subcategories. In order to make this research paper easy-to-follow for researchers at different levels (i.e., students, researchers, and senior researchers), we aimed at showing the evolution of each category over the decades rather than including many techniques with an exhaustive reference list. The experimental section was designed in a way to compare the performances of the techniques in two ways: (1) between all the categories (i.e., shallow unsupervised, shallow supervised, and deep feature extraction approaches) and (2) within each category by analyzing the corresponding subcategories. In this manner, a variety of subcategories is investigated detailing the evolution of shallow unsupervised approaches (i.e., conventional data projection schemes, band clustering/splitting techniques, low-rank reconstruction techniques, and manifold learning

techniques), shallow supervised approaches (i.e., class separation discriminant analysis, graph embedding discriminant analysis, regression-based representation learning, and joint & progressive learning strategy), and deep approaches (i.e., AE, CNN, RNN, and integrative approaches). Three recent hyperspectral datasets have been studied in this paper and the results are evaluated in terms of classification accuracies and the quality of the classification maps. The experiments carried out in this study showed that in terms of classification accuracies; 1) deep learning FE techniques (i.e., CNN and PCNN) can outperform the shallow ones in particular when a sufficient amount of training data is available 2) applying a dimensionality reduction step (such as PCA) prior to the deep learning techniques considerably improves their performances 3) shallow UFE techniques not only outperform the SFE ones but also are very competitive compared with deep FE ones. However, we should mention that the conclusions are limited by the experiments carried out on the three HSI datasets. In addition, this paper provides an impressive amount of codes and libraries mostly written in Python and Matlab to ease out the task of researchers in this very vibrant field of research.

# ACKNOWLEDGMENT

We would like to thank Prof. Melba Crawford for providing the Indian Pines 2010 Data and the National Center for Airborne Laser Mapping, the Hyperspectral Image Analysis Laboratory at the University of Houston, and the IEEE GRSS Image Analysis and Data Fusion Technical Committee. This work is partially supported by an Alexander von Humboldt research grant. We also would like to thank the AXA Research Fund for supporting the work of Prof. Jocelyn Chanussot. The corresponding author of this paper is Dr. Danfeng Hong.

# REFERENCES

[1] P. Ghamisi, J. Plaza, Y. Chen, J. Li, and A. J. Plaza, "Advanced spectral classifiers for hyperspectral images: A review," *IEEE Geosci. Remote Sens. Mag.*, vol. 5, no. 1, pp. 8-32, March 2017.

[2] D. W. Scott, *Multivariate Density Estimation*, John Wiley & Sons, New York, NY, 1992.

[3] E. J. Wegman, "Hyperdimensional data analysis using parallel coordinates," *J. Am. Stat. Assoc.*, vol. 85, no. 411, pp. 664-675, 1990.

[4] L. Jimenez and D. Landgrebe, "Supervised classification in high-dimensional space: geometrical, statistical, and asymptotical properties of multivariate data," *IEEE Trans. Syst. Man Cybern. C: Appl. Rev.*, vol. 28, no. 1, pp. 39-54, 1998.

[5] J. A. Benediktsson and P. Ghamisi, *Spectral-Spatial Classification of Hyperspectral Remote Sensing Images*. Artech House Publishers, INC, Boston, USA, 2015.

[6] K. Fukunaga, *Introduction to Statistical Pattern Recognition*, 2nd ed. Academic Press, Inc., San Diego, CA, 1990.

[7] D. A. Landgrebe, *Signal Theory Methods in Multispectral Remote Sensing*. John Wiley & Sons, Hoboken, NJ, 2003.

[8] G. Hughes, "On the mean accuracy of statistical pattern recognizers," *IEEE Trans. Inf. Theory*, vol. 14, pp. 55-63, 1968.

[9] Y. Qian, F. Yao, and S. Jia, "Band selection for hyperspectral imagery using affinity propagation," *IET Comput. Vis.*, vol. 3, no. 4, p. 213, 2009.

[10] L. Guanter and et. al., "The EnMAP spaceborne imaging spectroscopy mission for earth observation," *Remote Sens.*, vol. 7, no. 7, pp. 8830-8857, July 2015.
[11] P. Ghamisi, E. Maggiori, S. Li, R. Souza, Y. Tarablaka, G. Moser, A. De Giorgi, L. Fang, Y. Chen, M. Chi, S. B. Serpico, and J. A. Benediktsson, "New frontiers in spectral-spatial hyperspectral image classification: The latest advances based on mathematical morphology, markov random fields, segmentation, sparse representation, and deep learning," *IEEE Geosci. Remote Sens. Mag.*, vol. 6, no. 3, pp. 10–43, Sep. 2018.

[12] L. Zhang, Q. Zhang, B. Du, X. Huang, Y. Y. Tang, and D. Tao, "Simultaneous spectral-spatial feature selection and extraction for hyperspectral images," *IEEE Transactions on Cybernetics*, vol. 48, no. 1, pp. 16–28, Jan 2018.

[13] M. Fauvel, Y. Tarabalka, J. A. Benediktsson, J. Chanussot, and J. C. Tilton, "Advances in spectral-spatial classification of hyperspectral images," *Proceedings of the IEEE*, vol. 101, no. 3, pp. 652–675, 2013.

[14] W. Li, F. Feng, H. Li, and Q. Du, "Discriminant analysis-based dimension reduction for hyperspectral image classification: A survey of the most recent advances and an experimental comparison of different techniques," *IEEE Geosci. Remote Sens. Mag.*, vol. 6, no. 1, pp. 15–34, March 2018.

[15] X. Jia, B.-C. Kuo, and M. Crawford, "Feature mining for hyperspectral image classification," *Proc. IEEE*, vol. 101, no. 3, pp. 676–697, March 2013.

[16] W. Sun and Q. Du, "Hyperspectral band selection: A review," *IEEE Geosci. Remote Sens. Mag.*, vol. 7, no. 2, pp. 118–139, June 2019.

[17] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, and J. Chanussot, "Hyperspectral unmixing overview: Geometrical, statistical, and sparse regression-based approaches," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 5, no. 2, pp. 354–379, April 2012.

[18] P. Ghamisi, M. Dalla Mura, and J. A. Benediktsson, "A survey on spectralspatial classification techniques based on attribute profiles," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 5, pp. 2335–2353, May 2015.

[19] I. Jolliffe, *Principal Component Analysis*. Springer Verlag, 1986.

[20] J. Senthilnath, S. N. Omkar, V. Mani, N. Karnwal, and S. P. B., "Crop stage classification of hyperspectral data using unsupervised techniques," *IEEE J. Sel. Top. App. Earth Observ. Remote Sens.*, vol. 6, no. 2, pp. 861–866, April 2013.

[21] A. Green, M. Berman, P. Switzer, and M. Craig, "A transformation for ordering multispectral data in terms of image quality with implications for noise removal," *IEEE Trans. Geos. Remote Sens.*, vol. 26, no. 1, pp. 65–74, 1988.

[22] J. Lee, A. Woodyatt, and M. Berman, "Enhancement of high spectral resolution remote-sensing data by a noise-adjusted principal components transform," *IEEE Trans. Geos. Remote Sens.*, vol. 28, no. 3, pp. 295–304, 1990.

[23] A. Hyvärinen, J. Karhunen, and E. Oja, *Independent Component Analysis*, ser. Adaptive and Learning Systems for Signal Processing, Communications and Control Series. Wiley, 2001.

[24] A. Villa, J. Chanussot, C. Jutten, J. Benediktsson, and S. Moussaoui, "On the use of ICA for hyperspectral image analysis," in *Proc. IGARSS*, vol. 4, July 2009, pp. IV–97–IV–100.

[25] A. A. Nielsen, "Kernel maximum autocorrelation factor and minimum noise fraction transformations," *IEEE Trans. Image Process.*, vol. 20, no. 3, pp. 612–624, Mar 2011.

[26] F. Mei, C. Zhao, L. Wang, and H. Huo, "Anomaly detection in hyperspectral imagery based on kernel ICA feature extraction," in *Proceedings of the 2008 Second International Symposium on Intelligent Information Technology Application - Volume 1*, ser. IITA '08. Washington, DC, USA: IEEE Computer Society, 2008, pp. 869–873. [Online]. Available: http://dx.doi.org/10.1109/IITA.2008.98

[27] B. Schölkopf, A. J. Smola, and K.-R. Müller, "Advances in kernel methods," pp. 327–352, 1999.

[28] S. Marchesi and L. Bruzzone, "ICA and kernel ICA for change detection in multispectral remote sensing images," in *Proc. IGARSS*, vol. 2, July 2009, pp. II–980–II–983.

[29] M. Fauvel, J. Chanussot, and J. Benediktsson, "Kernel principal component analysis for the classification of hyperspectral remote sensing data over urban areas," *EURASIP J. Adv. Signal Process.*, vol. 2009, no. 1, p. 783194, 2009.

[30] L. M. Bruce and C. H. K. and, "Dimensionality reduction of hyperspectral data using discrete wavelet transform feature extraction," *IEEE Trans. Geosci. Remote Sens.*, vol. 40, no. 10, pp. 2331–2338, Oct 2002.

[31] S. Kumar, J. Ghosh, and M. M. Crawford, "Best-bases feature extraction algorithms for classification of hyperspectral data," *IEEE Trans. Geosci. Remote Sens.*, vol. 39, no. 7, pp. 1368–1379, July 2001.

[32] A. Martnez-UsMartinez-Uso, F. Pla, J. M. Sotoca, and P. Garca-Sevilla, "Clustering-based hyperspectral band selection using information measures," *IEEE Trans. Geosci. Remote Sens.*, vol. 45, no. 12, pp. 4158–4171, Dec 2007.

[33] C. Cariou, K. Chehdi, and S. Le Moan, "Bandclust: An unsupervised band reduction method for hyperspectral remote sensing," *IEEE Geosci. Remote Sens. Let.*, vol. 8, no. 3, pp. 565–569, May 2011.

[34] S. Rashwan and N. Dobigeon, "A split-and-merge approach for hyperspectral band selection," *IEEE Geosci. Remote Sens. Let.*, vol. 14, no. 8, pp. 1378–1382, Aug 2017.

[35] X. Kang, S. Li, and J. A. Benediktsson, "Feature extraction of hyperspectral images with image fusion and recursive filtering," *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 6, pp. 3742–3752, June 2014.

[36] X. Kang, S. Li, L. Fang, and J. A. Benediktsson, "Intrinsic image decomposition for feature extraction of hyperspectral images," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 4, pp. 2241–2253, April 2015.

[37] X. Jin, Y. Gu, and T. Liu, "Intrinsic image recovery from remote sensing hyperspectral images," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 1, pp. 224–238, Jan. 2019.

[38] X. Jin and Y. Gu, "Superpixel-based intrinsic image decomposition of hyperspectral images," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 8, pp. 4285–4295, Aug. 2017.

[39] P. Duan, X. Kang, S. Li, and P. Ghamisi, "Noise-robust hyperspectral image classification via multi-scale total variation," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 12, no. 6, pp. 1948–1962, 2019.

[40] B. Rasti, "Sparse Hyperspectral Image Modeling and Restoration," Ph.D. dissertation, Dep. of Electrical and Computer Engineering, University of Iceland, Dec. 2014.

[41] B. Rasti, J. Sveinsson, and M. Ulfarsson, "Wavelet-based sparse reduced-rank regression for hyperspectral image restoration," *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 10, pp. 6688–6698, Oct 2014.

[42] B. Rasti, M. O. Ulfarsson, and J. R. Sveinsson, "Hyperspectral feature extraction using total variation component analysis," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 12, pp. 6976–6985, Dec 2016.

[43] B. Rasti, P. Ghamisi, and M. O. Ulfarsson, "Hyperspectral feature extraction using sparse and smooth low-rank analysis," *Remote Sens.*, vol. 11, no. 2, p. 121, 2019. [Online]. Available: https://doi.org/10.3390/rs11020121

[44] J. B. Tenenbaum, V. de Silva, and J. C. Langford, "A global geometric framework for nonlinear dimensionality reduction," *Science*, vol. 290, no. 5500, p. 2319, 2000.

[45] C. M. Bachmann, T. L. Ainsworth, and R. A. Fusina, "Exploiting manifold geometry in hyperspectral imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 43, no. 3, pp. 441–454, March 2005.

[46] S. T. Roweis and L. K. Saul, "Nonlinear dimensionality reduction by locally linear embedding," *Science*, vol. 290, pp. 2323–2326, 2000.

[47] M. Belkin and P. Niyogi, "Laplacian eigenmaps and spectral techniques for embedding and clustering," in *Proc. NIPS*. MIT Press, 2002, pp. 585–591.

[48] X. He and P. Niyogi, "Locality Preserving Projections," in *Proc. NIPS*, S. Thrun, L. Saul, and B. Scholkopf, Eds. MIT Press, 2003.

[49] S. Yan, D. Xu, B. Zhang, H.-J. Zhang, Q. Yang, and S. Lin, "Graph embedding and extensions: A general framework for dimensionality reduction," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 29, no. 1, pp. 40–51, Jan 2007.

[50] X. Wu, D. Hong, J. Tian, J. Chanussot, W. Li, and R. Tao, "Orsim detector: A novel object detection framework in optical remote sensing imagery using spatial-frequency channel features," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 7, pp. 5146–5158, 2019.

[51] R. O. Duda and P. E. Hart, *Pattern classification and scene analysis*. Wiley New York, 1973, vol. 3.

[52] R. Hang, Q. Liu, H. Song, and Y. Sun, "Matrix-based discriminant subspace ensemble for hyperspectral image spatial-spectral feature fusion," *IEEE Trans. Geos. Remote Sens.*, vol. 54, no. 2, pp. 783–794, 2016.

[53] C. Lee and D. Landgrebe, "Feature extraction based on decision boundaries," *IEEE Trans. Pattern Anal. Mach. Intell.*, vol. 15, no. 4, pp. 388–400, 1993.

[54] B. Guo, S. R. Gunn, R. I. Damper, and J. D. Nelson, "Band selection for hyperspectral image classification using mutual information," *IEEE Geosci. Remote Sens. Let.*, vol. 3, no. 4, pp. 522–526, 2006.

[55] H. Shi, Y. Shen, and Z. Liu, "Hyperspectral bands reduction based on rough sets and fuzzy c-means clustering," in *IEEE Instrumentation Technology Conference (IMTC)*, vol. 2. IEEE, 2003, pp. 1053–1056.

[56] J. Yang and J. Yang, "Why can Ida be performed in pca transformed space?" *Pattern Recog.*, vol. 36, no. 2, pp. 563–566, 2003.
[57] T. V. Bandos, L. Bruzzone, and G. Camps-Valls, "Classification of hyperspectral images with regularized linear discriminant analysis," *IEEE Trans. Geos. Remote Sens.*, vol. 47, no. 3, pp. 862–873, 2009.

[58] M. Sugiyama, "Dimensionality reduction of multimodal labeled data by local fisher discriminant analysis," *J. Mach. Learn. Res.*, vol. 8, no. May, pp. 1027–1061, 2007.

[59] M. Imani and H. Ghasseman, "Feature space discriminant analysis for hyperspectral data feature reduction," *ISPRS J. Photogramm. Remote Sens.*, vol. 102, pp. 1–13, 2015.

[60] S. Patra, P. Modi, and L. Bruzzone, "Hyperspectral band selection based on rough set," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 10, pp. 5495–5503, 2015.

[61] X. Cao, T. Xiong, and L. Jiao, "Supervised band selection using local spatial information for hyperspectral image," *IEEE Geosci. Remote Sens. Lett.*, vol. 13, no. 3, pp. 329–333, 2016.

[62] M. Belkin and P. Niyogi, "Laplacian eigenmaps for dimensionality reduction and data representation," *Neural Comput.*, vol. 15, no. 6, pp. 1373–1396, 2003.

[63] X. He and P. Niyogi, "Locality preserving projections," in *Proc. NIPS*, 2004, pp. 153–160.

[64] L. Zelnik-Manor and P. Perona, "Self-tuning spectral clustering," in *Proc. NIPS*, 2005, pp. 1601–1608.

[65] N. H. Ly, Q. Du, and J. E. Fowler, "Sparse graph-based discriminant analysis for hyperspectral imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 52, no. 7, pp. 3872–3884, 2014.

[66] N. H. Ly, Q. Du, and J. E. Fowler, "Collaborative graph-based discriminant analysis for hyperspectral imagery," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 7, no. 6, pp. 2688–2696, 2014.

[67] D. Hong, N. Yokoya, and X. X. Zhu, "Learning a robust local manifold representation for hyperspectral dimensionality reduction," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 10, no. 6, pp. 2960–2975, 2017.

[68] D. Hong, N. Yokoya, J. Chanussot, J. Xu, and X. X. Zhu, "Learning to propagate labels on graphs: An iterative multitask regression framework for semi-supervised hyperspectral dimensionality reduction," *ISPRS J. Photogramm. Remote Sens.*, vol. 158, pp. 35–49, 2019.

[69] L. Zhang, M. Yang, and X. Feng, "Sparse representation or collaborative representation: Which helps face recognition?" in *Proc. ICCV*. IEEE, 2011, pp. 471–478.

[70] D. Hong and X. X. Zhu, "SULoRA: Subspace unmixing with low-rank attribute embedding for hyperspectral data analysis," *IEEE J. Sel. Topics Signal Process.*, vol. 12, no. 6, pp. 1351–1363, 2018.

[71] D. Hong, N. Yokoya, J. Chanussot, and X. X. Zhu, "An augmented linear mixing model to address spectral variability for hyperspectral unmixing," *IEEE Trans. Image Process.*, vol. 28, no. 4, pp. 1923–1938, 2019.

[72] H. Huang, F. Luo, J. Liu, and Y. Yang, "Dimensionality reduction of hyperspectral images based on sparse discriminant manifold embedding," *ISPRS J. Photogramm. Remote Sens.*, vol. 106, pp. 42–54, 2015.

[73] Z. Xue, P. Du, J. Li, and H. Su, "Simultaneous sparse graph embedding for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 53, no. 11, pp. 6114–6133, 2015.

[74] W. Li and Q. Du, "Laplacian regularized collaborative graph for discriminant analysis of hyperspectral imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 12, pp. 7066–7076, 2016.

[75] W. Li, J. Liu, and Q. Du, "Sparse and low-rank graph for discriminant analysis of hyperspectral imagery," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 7, pp. 4094–4105, 2016.

[76] L. Pan, H.-C. Li, Y.-J. Deng, F. Zhang, X.-D. Chen, and Q. Du, "Hyperspectral dimensionality reduction by tensor sparse and low-rank graph-based discriminant analysis," *Remote Sens.*, vol. 9, no. 5, p. 452, 2017.

[77] K.-R. Müller, S. Mika, G. Rätsch, K. Tsuda, and B. Schölkopf, "An introduction to kernel-based learning algorithms," *IEEE Trans. Neural Netw.*, vol. 12, no. 2, pp. 1–14.

[78] S. Gao, I. W.-H. Tsang, and L.-T. Chia, "Kernel sparse representation for image classification and face recognition," in *Proc. ECCV*. Springer, 2010, pp. 1–14.

[79] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, "Distributed optimization and statistical learning via the alternating direction method of multipliers," *Found. Trends Mach. Learn.*, vol. 3, no. 1, p. 1122, Jan. 2011. [Online]. Available: https://doi.org/10.1561/2200000016

[80] G. H. Golub, P. C. Hansen, and D. P. O'Leary, "Tikhonov regularization and total least squares," *SIAM Journal on Matrix Analysis and Applications*, vol. 21, no. 1, pp. 185–194, 1999. [Online]. Available: https://doi.org/10.1137/S0895479897326432

[81] K.-C. Li, "Sliced inverse regression for dimension reduction," *J. Am. Stat. Assoc.*, vol. 86, no. 414, pp. 316–327, 1991.

[82] T. Suzuki and M. Sugiyama, "Sufficient dimension reduction via squared-loss mutual information estimation," *Neural Comput.*, vol. 25, no. 3, pp. 725–758, 2013.

[83] J. Sainui and M. Sugiyama, "Direct approximation of quadratic mutual information and its application to dependence-maximization clustering," *IEICE Trans. Inf. Syst.*, vol. 96, no. 10, pp. 2282–2285, 2013.

[84] V. Tangkaratt, H. Sasaki, and M. Sugiyama, "Direct estimation of the derivative of quadratic mutual information with application in supervised dimension reduction," *Neural Comput.*, vol. 29, no. 8, pp. 2076–2122, 2017.

[85] D. Hong, N. Yokoya, J. Chanussot, and X. X. Zhu, "Cospace: Common subspace learning from hyperspectral-multispectral correspondences," *IEEE Trans. Geos. Remote Sens.*, vol. 57, no. 7, pp. 4349–4359, 2019.

[86] D. Hong, N. Yokoya, N. Ge, J. Chanussot, and X. X. Zhu, "Learnable manifold alignment (LeMA): A semi-supervised cross-modality learning framework for land cover and land use classification," *ISPRS J. Photogramm. Remote Sens.*, vol. 147, pp. 193–205, 2019.

[87] D. Hong, N. Yokoya, J. Xu, and X. Zhu, "Joint & progressive learning from high-dimensional data for multi-label classification," in *Proc. ECCV*, 2018, pp. 469–484.

[88] L. Zhang, L. Zhang, and B. Du, "Deep learning for remote sensing data: A technical tutorial on the state of the art," *IEEE Geosci. Remote Sens. Mag.*, vol. 4, no. 2, pp. 22–40, 2016.

[89] S. Li, W. Song, L. Fang, Y. Chen, P. Ghamisi, and J. A. Benediktsson, "Deep learning for hyperspectral image classification: An overview," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 9, pp. 6690–6709, Sep. 2019.

[90] G. Licciardi and J. Chanussot, "Spectral transformation based on nonlinear principal component analysis for dimensionality reduction of hyperspectral images," *Eur. J. Remote Sens.*, vol. 51, no. 1, pp. 375–390, 2018.

[91] Y. Chen, Z. Lin, X. Zhao, G. Wang, and Y. Gu, "Deep learning-based classification of hyperspectral data," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 7, no. 6, pp. 2094–2107, 2014.

[92] C. Tao, H. Pan, Y. Li, and Z. Zou, "Unsupervised spectral-spatial feature learning with stacked sparse autoencoder for hyperspectral imagery classification," *IEEE Geosci. Remote Sens. Lett.*, vol. 12, no. 12, pp. 2438–2442, 2015.

[93] R. Kemker and C. Kanan, "Self-taught feature learning for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 5, pp. 2693–2705, 2017.

[94] X. Kang, C. Li, S. Li, and H. Lin, "Classification of hyperspectral images by gabor filtering based deep network," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 4, pp. 1166–1178, 2018.

[95] C. Deng, Y. Xue, X. Liu, C. Li, and D. Tao, "Active transfer learning network: A unified deep joint spectral-spatial feature learning model for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 3, pp. 1741–1754, 2019.

[96] X. Ma, H. Wang, and J. Geng, "Spectral-spatial classification of hyperspectral image based on deep auto-encoder," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 9, no. 9, pp. 4073–4085, 2016.

[97] X. Sun, F. Zhou, J. Dong, F. Gao, Q. Mu, and X. Wang, "Encoding spectral and spatial context information for hyperspectral image classification," *IEEE Geosci. Remote Sens. Lett.*, vol. 14, no. 12, pp. 2250–2254, 2017.

[98] P. Zhou, J. Han, G. Cheng, and B. Zhang, "Learning compact and discriminative stacked autoencoder for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, 2019.

[99] W. Hu, Y. Huang, L. Wei, F. Zhang, and H. Li, "Deep convolutional neural networks for hyperspectral image classification," *Sensors*, vol. 2015, 2015.

[100] W. Zhao and S. Du, "Spectral-spatial feature extraction for hyperspectral image classification: A dimension reduction and deep learning approach," *IEEE Trans. Geosci. Remote Sens.*, vol. 54, no. 8, pp. 4544–4554, 2016.

[101] 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.

[102] X. Liu, C. Deng, J. Chanussot, D. Hong, and B. Zhao, "Stfnet: A two-stream convolutional neural network for spatiotemporal image fusion," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 9, pp. 6552–6564, 2019.

[103] W. Song, S. Li, L. Fang, and T. Lu, "Hyperspectral image classification with deep feature fusion network," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 6, pp. 3173–3184, 2018.
[104] X. Ma, A. Fu, J. Wang, H. Wang, and B. Yin, "Hyperspectral image classification based on deep deconvolution network with skip architecture," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 8, pp. 4781-4791, 2018.

[105] P. Ghamisi, Y. Chen, and X. X. Zhu, "A self-improving convolution neural network for the classification of hyperspectral data," *IEEE Geosci. Remote Sens. Lett.*, vol. 13, no. 10, pp. 1537-1541, Oct 2016.

[106] Y. Zhan, D. Hu, H. Xing, and X. Yu, "Hyperspectral band selection based on deep convolutional neural network and distance density," *IEEE Geosci. Remote Sens. Lett.*, vol. 14, no. 12, pp. 2365-2369, Dec 2017.

[107] Y. Kong, X. Wang, and Y. Cheng, "Spectral-spatial feature extraction for hsi classification based on supervised hypergraph and sample expanded cnn," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 11, pp. 4128-4140, Nov 2018.

[108] J. Yang, Y. Zhao, and J. C. Chan, "Learning and transferring deep joint spectralspatial features for hyperspectral classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 8, pp. 4729-4742, Aug 2017.

[109] S. Mei, J. Ji, J. Hou, X. Li, and Q. Du, "Learning sensor-specific spatial-spectral features of hyperspectral images via convolutional neural networks," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 8, pp. 4520-4533, Aug 2017.

[110] L. Jiao, M. Liang, H. Chen, S. Yang, H. Liu, and X. Cao, "Deep fully convolutional network-based spatial distribution prediction for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 10, pp. 5585-5599, Oct 2017.

[111] G. Cheng, Z. Li, J. Han, X. Yao, and L. Guo, "Exploring hierarchical convolutional features for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 11, pp. 6712-6722, Nov 2018.

[112] M. Liang, L. Jiao, S. Yang, F. Liu, B. Hou, and H. Chen, "Deep multiscale spectral-spatial feature fusion for hyperspectral images classification," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, vol. 11, no. 8, pp. 2911-2924, Aug 2018.

[113] H. Wu and S. Prasad, "Semi-supervised deep learning using pseudo labels for hyperspectral image classification," *IEEE Trans. Image Process.*, vol. 27, no. 3, pp. 1259-1270, 2018.

[114] X. Xu, W. Li, Q. Ran, Q. Du, L. Gao, and B. Zhang, "Multisource remote sensing data classification based on convolutional neural network," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 2, pp. 937-949, 2018.

[115] Z. Zhong, J. Li, Z. Luo, and M. Chapman, "Spectral-spatial residual network for hyperspectral image classification: A 3-d deep learning framework," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 2, pp. 847-858, 2018.

[116] Q. Liu, F. Zhou, R. Hang, and X. Yuan, "Bidirectional-convolutional lstm based spectral-spatial feature learning for hyperspectral image classification," *Remote Sens.*, vol. 9, no. 12, p. 1330, 2017.

[117] F. Zhou, R. Hang, Q. Liu, and X. Yuan, "Hyperspectral image classification using spectral-spatial lstms," *Neurocomputing*, vol. 328, pp. 39-47, 2019.

[118] X. Zhang, Y. Sun, K. Jiang, C. Li, L. Jiao, and H. Zhou, "Spatial sequential recurrent neural network for hyperspectral image classification," *IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens.*, no. 99, pp. 1-15, 2018.

[119] C. Shi and C.-M. Pun, "Multi-scale hierarchical recurrent neural networks for hyperspectral image classification," *Neurocomputing*, vol. 294, pp. 82-93, 2018.

[120] A. Ma, A. M. Filippi, Z. Wang, and Z. Yin, "Hyperspectral image classification using similarity measurements-based deep recurrent neural networks," *Remote Sens.*, vol. 11, no. 2, p. 194, 2019.

[121] Y. Xu, L. Zhang, B. Du, and F. Zhang, "Spectral-spatial unified networks for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 10, pp. 5893-5909, 2018.

[122] R. Hang, Q. Liu, D. Hong, and P. Ghamisi, "Cascaded recurrent neural networks for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 57, no. 8, pp. 5384-5394, 2019.

[123] S. Hao, W. Wang, Y. Ye, T. Nie, and L. Bruzzone, "Two-stream deep architecture for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 56, no. 4, pp. 2349-2361, 2018.

[124] P. Vincent, H. Larochelle, I. Lajoie, Y. Bengio, and P. Manzagol, "Stacked denoising autoencoders," *J. Mach. Learn. Res.*, vol. 11, no. 12, pp. 3371-3408, 2010.

[125] A. Krizhevsky, I. Sutskever, and G. Hinton, "Imagenet classification with deep convolutional neural networks," in *Proc. NIPS*, Lake Tahoe, Nevada, USA, 2012, pp. 1527-1554.

[126] R. Kemker and C. Kanan, "Self-taught feature learning for hyperspectral image classification," *IEEE Trans. Geosci. Remote Sens.*, vol. 55, no. 5, pp. 2693-2705, 2017.

[127] P. Ghamisi, N. Yokoya, J. Li, W. Liao, S. Liu, J. Plaza, B. Rasti, and A. Plaza, "Advances in hyperspectral image and signal processing: A comprehensive overview of the state of the art," *IEEE Geosci. Remote Sens. Mag.*, vol. 5, no. 4, pp. 37-78, Dec 2017.

[128] B. Rasti, M. Ulfarsson, and J. Sveinsson, "Hyperspectral subspace identification using SURE," *IEEE Geosci. Remote Sens. Let.*, vol. 12, no. 12, pp. 2481-2485, Dec. 2015.

[129] J. Bioucas-Dias and J. Nascimento, "Hyperspectral subspace identification," *IEEE Trans. Geos. Remote Sens.*, vol. 46, no. 8, pp. 2435-2445, 2008.

[130] B. Rasti, P. Scheunders, P. Ghamisi, G. Licciardi, and J. Chanussot, "Noise reduction in hyperspectral imagery: Overview and application," *Remote Sens.*, vol. 10, no. 3, p. 482, 2018. [Online]. Available: http://www.mdpi.com/2072-4292/10/3/482