Hyperspectral Remote Sensing Image Classification Based on Partitioned Random Projection Algorithm

Abstract

Dimensionality reduction based on random projection (RP) includes two problems, namely, the dimensionality is limited by the data size and the class separability of the dimensionality reduction results is unstable due to the randomly generated projection matrix. These problems make the RP algorithm unsuitable for large-size hyperspectral image (HSI) classification. To solve these problems, this paper presents a new partitioned RP (PRP) algorithm and proves its rationality in theory. First, a large-size HSI is evenly divided into multiple small-size sub-HSIs. Afterwards, the projection matrix that maximizes the class separability is selected from multiple samplings in which the class dissimilarity measurement is defined as large inter-class distance and small intra-class variance. By using the same projection matrix, each small-size sub-HSI is projected to generate a low dimensional sub-HSI, thereby generating a low dimensional HSI. Next, the minimum distance (MD) classifier is utilized to classify the low dimensional HSI obtained by the PRP algorithm. Finally, four real HSIs are used for experiments, and three of the most popular classification algorithms based on RP are selected as comparison algorithms to validate the effectiveness of the proposed algorithm. The classification performance is evaluated with the kappa coefficient, overall accuracy (OA), average accuracy (AA), average precision rate (APR), and running time. Experimental results indicate that the proposed algorithm can obtain reliable classification results in a very short time.

1. Introduction

The wavelength range of hyperspectral images (HSIs) is continuous and dense from visible light to short-wave infrared [1]. Compared with panchromatic images and multispectral remote sensing images, HSIs contain richer spectral information and finer texture information of ground objects, which provide a solid data foundation for high-precision ground object classification [2,3]. Therefore, more and more researchers apply it to the fields of medicine [4,5] and agriculture [6,7]. However, for commonly configured hardware, the classification of large-size HSIs has some problems, such as long running time and a large amount of data. The most straightforward solution is to perform dimensionality reduction before HSI classification, thereby reducing the computational load through dimensionality reduction [8,9].

Popular HSI dimensionality reduction methods are roughly divided into two types: band selection [10,11] and feature extraction [12,13]. The band selection method directly discards most of the bands and only selects a part of the bands for subsequent HSI analysis. Currently, Zhang et al. [14] proposed a marginalized graph self-representation (MGSR) method for unsupervised hyperspectral band selection, which considers the differences between different homogenous regions. First, super-pixel segmentation is used to construct the structure map of HSI [15]. Additionally, considering the relationship between adjacent pixels of the same segmentation, the damage is marginalized by an alternating optimization algorithm to search for the best band selection scheme. Computation of the graph matrix takes a lot of time, and the result of super-pixel segmentation affects the result of subsequent band selection. Mou et al. [16] studied deep reinforcement learning models for HSI analysis and first proposed an automatic learning strategy for band selection under the framework of the Markov decision. The method first used a parametric method to rank the bands and then utilized a deep reinforcement learning model to select the best subset of bands. Because it is a learning-based algorithm, it takes much time during the training phase to explore efficient band selection strategies. Feature extraction generally uses mapping to linearly transform all bands of HSIs. Tamilarasi et al. [17] proposed a new dimensionality reduction and classification technique, which extracted and classified road and building features from HSIs with high accuracy. This technique combined independent component analysis (ICA) [18], principal component analysis (PCA) [19], fully convolutional network (FCN) [20], and support vector machine (SVM) models [21]. The above dimensionality reduction method is cumbersome and cannot be calculated quickly.

The random projection (RP) algorithm is an emerging dimensionality reduction technique that satisfies the Johnson-Lindenstrauss (JL) lemma [22,23,24]. The algorithm does not require any prior knowledge of the original data [25]. In addition, a randomly generated projection matrix is used to achieve dimensionality reduction, where the column vector of the projection matrix is unit length. This makes the RP algorithm easier to implement and more convenient to calculate [26]. Additionally, the RP algorithm can guarantee that the distance between any pair of vectors before and after projection is kept within a small range [27,28]. Because many machine learning algorithms use the distance information between vectors to operate and execute [29,30], the dataset obtained after dimensionality reduction by the RP algorithm can be well applied to machine learning algorithms. Therefore, more and more researchers apply it to the dimensionality reduction of HSI data.

At present, many researchers interested in the field of the RP algorithm have proposed many classification algorithms based on the RP algorithm [31,32]. These methods can be grouped into three types: combined feature extraction method [33], separability-boosting method [34], and ensemble method [35]. Combined feature extraction methods combine RP and other feature extraction methods. Zhao and Mao [36] proposed a semi-random projection method, which uses linear discriminant analysis (LDA) to calculate each column vector in the projection matrix and calculates the projection matrix by repeating it multiple times. It strikes a good balance between computational complexity and classification accuracy. This method can effectively preserve the distance invariance of the data when dealing with high-dimensional data. Deegalla and Bostrom [37] compared the dimensionality reduction effect of PCA and the RP algorithm and combined PCA and RP to propose a nearest-neighbor classification algorithm. This method not only improves efficiency but also enables more efficient nearest-neighbor classification. However, the method cannot achieve dimensionality reduction when dealing with large-size HSIs. Separability-boosting methods preserve the class separability of the original data during projection. Zhao et al. [38] proposed a new tighter RP with minimal intra-class variance (TRP-MIV) algorithm for HSI classification, which supplied a promising way for HSI dimensionality reduction. Nevertheless, when the number of hyperspectral vectors is very large, the dimensionality reduction effectiveness of the TRP-MIV algorithm is still limited. The ensemble methods combine RP and the fuzzy clustering algorithm. Rathore et al. [39] studied the cumulative agreement fuzzy c-means (CAFCM) algorithm, which uses RP for generating multiple dimensionality reduction results. In addition, Anderlucci et al. [40] proposed a model-based high dimensional data clustering method based on the idea of ensemble methods. The method first generates a set of low dimensional independent projection results and performs model-based clustering on each of them. Then, the projection of the best grouping structure is displayed, and the final partition is obtained by aggregating the clusters found in the projection by consensus. However, the dimensionality reduction results produced by each projection are only guaranteed to be like the original data and cannot be applied to classification tasks. The RP algorithm has the following shortcomings in large-size HSI classification. (1) The lowest projection dimensionality of the RP algorithm is independent of the high dimensionality, which is positively correlated with the number of hyperspectral vectors [41]. The larger number of hyperspectral vectors, the larger the lowest projection dimensionality. (2) Different dimensionality reduction results will be generated by different projection matrices, which lead to the class separability of dimensionality reduction results being very unstable. Thus, the RP algorithm cannot play a dimensionality reduction effect on the HSI with a large number of hyperspectral vectors, which is not conducive to large-size HSI classification.

Aiming at the above problems existing in the RP algorithm in large-scale dimensionality reduction and classification, an HSI classification algorithm based on partitioned RP (PRP) is proposed. Firstly, the PRP algorithm is proposed by means of image division, and the distance preservation of the hyperspectral vector pair after partitioning is still satisfied from two perspectives within each sub-HSI and between any two different sub-HSIs. In addition, a classifier based on the PRP algorithm is designed, in which the optimal projection matrix is determined according to the principle of large inter-class distance and small intra-class variance, and minimum distance (MD) [42,43] is used to classify the low dimensional HSI received by the PRP algorithm. The experimental results show that the PRP algorithm can project the large-size HISs into a subspace with a lower dimensionality than the RP algorithm, and the proposed classification algorithm can obtain good classification results. The structure of this paper is organized as below. Section 2 and Section 3 give the materials and the proposed algorithm, respectively. The experimental results and discussion are provided in Section 4 and Section 5. Finally, this paper is concluded in Section 6.

2. Materials

This paper uses four publicly available datasets with validation data, namely, real Pavia Centre, Salinas, Chikusei, and LongKou [44,45] images. The four images are respectively located in Pavia Centre in northern Italy, Salinas in Northwestern Indiana, Chikusei in Japan, and LongKou in China. Details of these images are listed in

Table 1. Details of four real HSIs.

	Pavia Centre	Salinas	Chikusei	LongKou
URL for data source	http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes (accessed on 12 July 2021)	http://www.ehu.eus/ccwintco/index.php/Hyperspectral_Remote_Sensing_Scenes (accessed on 12 July 2021)	https://naotoyokoya.com/Download.html (accessed on 11 April 2016)	http://rsidea.whu.edu.cn/resource_WHUHi_sharing.htm (accessed on 1 October 2020)
Year	2003	1998	2014	2018
Sensor	ROSIS	AVIRIS	Headwall Hyperspec-VNIR-C	Nano-Hyperspec
Spatial resolution	1.3 m	3.7 m	2.5 m	0.463 m
Number of bands	102	204	128	270

, including download address, year, sensor, spatial resolution, and number of bands. Spatial resolution represents the actual distance on the ground represented by a pixel.

It is worth noting that this paper uses HSIs of sizes 1096 × 531, 512 × 127, 1000 × 800, and 550 × 400, with homogeneous regions of 9, 13, 7, and 9, respectively, as shown in Figure 1. The classes included in the Pavia Centre image are water, trees, meadows, self-blocking bricks, bare soil, asphalt, bitumen, tiles, and shadows. The classes included in the Salinas image are broccoli green weeds2, fallow, fallow rough plow, fallow smooth, stubble, celery, grapes untrained, soil vineyard develop, corn senesced weeds, lettuce romaine 4 week, lettuce romaine 5 week, lettuce romaine 6 week, and lettuce romaine 7 week. The classes included in the Chikusei image are water, bare soil (school), bare soil (farmland), natural plants, glass, rice field (grown), and row crops. The classes included in the LongKou image are corn, cotton, sesame, broad leaf soybean, narrow leaf soybean, rice, water, roads and houses, and mixed weed. Specifically, Figure 1 shows the false-color images, the standard images (i.e., validation data), and the legends. The band combinations used to replace the RGB bands in the false-color images of the four images are (49, 31, 15), (29, 20, 12), (54, 35, 22) and (130, 65, 18), respectively. Moreover, the specific classes are explained in the legends. In addition, because many vectors in these images do not contain any valid information, these vectors are discarded before applying. Therefore, the numbers of spectral vectors in these scenes are 109,794, 20,655, 9435, and 204,542, respectively.

3. The Proposed Algorithm

3.1. Random Projection

RP is an effective dimensionality reduction tool that can keep the distance between vector pairs before and after projection unchanged. The usual version of the RP algorithm is defined below [26].

Theorem 1.

Supposing that A is any S × D matrix, where S is the data size and D is the feature dimensionality. Its row vector is a_s in D dimensional spaces, where s is the index of row vectors and s = 1, …, S. For constants ε and β > 0, an integer K_RP is chosen as follows,

$$K_{\mathrm{RP}}\ge K_{\mathrm{RP}}^0=\frac{4+2\beta }{\frac{{\epsilon }^2}{2}-\frac{{\epsilon }^3}{3}}\log S$$

(1)

where K⁰_RP is the lowest projection dimensionality of the RP algorithm. R_RP is any D × K_RP matrix whose entry belongs to the standard normal Gaussian distribution. The S × K_RP matrix

$$\mathit{B}=\frac{1}{\sqrt{K_{\mathrm{RP}}}}\mathit{A}{\mathit{R}}_{\mathrm{RP}}$$

(2)

is the projection of A in K_RP dimensional subspaces. The row vector of matrix B is b_s. Then, with the probability of no distortion at least

P_RP = 1 − S ^−β

(3)

for any two row vectors a_s and a_s′ of matrix A, the distance of the row vectors b_s and b_s′ of matrix B is preserved as follows,

(1 − ε) ||a_s − a_s′||² ≤ ||b_s − b_s′||² ≤ (1 + ε) ||a_s − a_s′||²

(4)

where s′ is the index of row vectors.

It can be seen from the above theorem that RP is a computationally simple method. The order of forming a projection matrix R_RP and projecting the high dimensional data A into the K_RP dimensional subspaces is O (D × K_RP × S). If the number of vectors S becomes larger, the lowest projection dimensionality K⁰_RP will become larger according to Equation (1) and the probability of no distortion P_RP will also become larger according to Equation (3). Therefore, when the number of vectors is large, although the probability of no distortion is high, the low dimensionality K_RP may be higher than the high dimensionality D, so the dimensionality reduction effect cannot be achieved. Besides, the generation of the projection matrix is completely random, which will lead to instability in the class separability of the generated dimensionality reduction result. In simpler terms, it cannot guarantee that the class separability of the dimensionality reduction result is good for the classification task.

3.2. Partitioned Random Projection

The lowest projection dimensionality K⁰_RP of the RP algorithm is strictly limited by the parameters ε, β, and the number of hyperspectral vectors S. The Pavia Centre image with 102 bands is taken as an example, of which the number of hyperspectral vectors is 109,794. When the fixed parameters ε = 1 and β = 0.5, according to Equation (1), the lowest projection dimensionality of the RP algorithm is 349, which is much higher than the number of bands, which is 102. Therefore, when processing large-size HSIs, if the number of bands is very low, the RP algorithm will not be able to achieve the dimensionality reduction of large-size HSIs. Considering this problem, a new PRP algorithm is introduced in this section. PRP, a dimensionality reduction algorithm, can project a large-size HSI into a low dimensional subspace without serious distortion on pairwise distance.

The HSI can be expressed as a matrix U = [u₁; …; u_s; …; u_S], where S is the number of hyperspectral vectors, s is the index of hyperspectral vectors, and u_s is the sth hyperspectral vector. The implementation of the PRP algorithm involves partitioning an HSI into multiple sub-HSIs evenly according to the horizontal direction and then projecting each sub-HSI into a low dimensional subspace. Specifically, the matrix of the HSI can also be defined as U = [U¹; …; U^m; …; U^M], where M is the number of sub-matrices, and m is the index of sub-matrices. Moreover, U^m = [u^m₁; …; u^m_n; …; u^m_N] is the mth N × D sub-matrix, where n is the index of hyperspectral vectors in the sub-matrix, and N is the number of hyperspectral vectors in the sub-matrix that is the same for all sub-matrices. Thus, S = M × N. Note that the dimensionality reduction method in this paper requires that an HSI must be divided equally. Otherwise, the first few hyperspectral vectors can be regarded as edges and discarded so that the HSI can be divided equally.

Next, the rationality of the distance preservation property guaranteed by this partition method is proven from two perspectives: within each sub-matrix and between any two different sub-matrices.

The distance changes before and after the projection of any two hyperspectral vectors within each sub-matrix are explained as follows. The conclusion of this theorem comes from the RP algorithm (Theorem 1).

Theorem 2.

For constants ε and β > 0, an integer K_PRP such that

$$K_{\mathrm{PRP}}\ge K_{\mathrm{PRP}}^0=\frac{4+2\beta }{\frac{{\epsilon }^2}{2}-\frac{{\epsilon }^3}{3}}\log N=\frac{4+2\beta }{\frac{{\epsilon }^2}{2}-\frac{{\epsilon }^3}{3}}\log \frac{S}{M}$$

(5)

where K⁰_PRP is the lowest projection dimensionality of the proposed algorithm. R_PRP is a D × K_PRP matrix whose entry belongs to the standard normal Gaussian distribution. The N × K_PRP sub-matrix

$${\mathit{V}}^m=\frac{1}{\sqrt{K_{\mathrm{PRP}}}}{\mathit{U}}^m{\mathit{R}}_{\mathrm{PRP}}$$

(6)

is the projection of sub-matrix U^m in K_PRP dimensional subspaces. The low dimensional hyperspectral vector of sub-matrix V^m is v^m_n. Namely, V^m = {v^m₁; …; v^m_n; …; v^m_N}. Then, according to Equation (3), with the probability of no distortion at least

P_PRP1 = 1 − N ^−β = 1 − (S/M)^−β

(7)

for any two hyperspectral vectors u^m_n and u^m_n′ of sub-matrix U^m, the distance of the two low dimensional hyperspectral vectors v^m_n and v^m_n′ of sub-matrix V^m is preserved based on Equation (4),

(1 − ε) ||u^m_n − u^m_n′||² ≤ ||v^m_n − v^m_n′||² ≤ (1 + ε) ||u^m_n − u^m_n′||²

(8)

where n′ is the index of hyperspectral vectors in the sub-matrix.

According to Equation (5), the lowest projection dimensionality K⁰_PRP decreases with the increase in the number of sub-matrices M and the decrease of parameter β. Moreover, with the increase of parameter ε, the change trend of lowest projection dimensionality K⁰_PRP is to first decrease and then increase. Comparing Equation (1) and Equation (5), the lowest projection dimensionality of the PRP algorithm is smaller than that of the RP algorithm under the same parameters. Note that when M is equal to one, it is the projection without partitioning, that is, the RP algorithm. When M is equal to S, there is only one hyperspectral vector in each sub-matrix. Meanwhile, the lowest projection dimensionality is zero, according to Equation (5), so it can be projected into any space with a dimensionality greater than zero.

A projection matrix is randomly generated, and each sub-matrix is projected into the same low dimensional space with this projection matrix. On this basis, the distance changes before and after the projection of any two hyperspectral vectors in any two different sub-matrices are introduced. It is described in detail as Theorem 3. (Proof of Theorem 3 can be found in Appendix A).

Theorem 3.

For constants ε and β > 0, an integer K_PRP is chosen according to Equation (5). R_PRP is a D × K_PRP matrix whose entries belong to the standard normal Gaussian distribution. The sub-matrix V^m is obtained based on Equation (6). Then, for the same ε, β, K_PRP, and R_PRP, with the probability of no distortion at least

P_PRP2 = 1 − 2N ^−β = 1 − 2 (S/M)^−β

(9)

for any hyperspectral vector u^m_n in sub-matrix U^m and hyperspectral vector u^m′_n′ in another sub-matrix U^m′, the distance of the low dimensional hyperspectral vector v^m_n in sub-matrix V^m and low dimensional hyperspectral vector v^m′_n′ in another sub-matrix V^m′ is preserved,

(1 − ε) ||u^m_n − u^m′_n′||² ≤ ||v^m_n − v^m′_n′||² ≤ (1 + ε) ||u^m_n − u^m′_n′||²

(10)

where m′ is the index of sub-matrices. Significantly, m′ is not equal to m (m′ ≠ m).

Combining Theorem 2 and Theorem 3, the distance before and after the projection of any two hyperspectral vectors in HSIs can also remain unchanged with an acceptable probability of no distortion. Specifically, the minimum value of probabilities of no distortion within each sub-matrix and between two different sub-matrices is taken as the probability of no distortion of the PRP algorithm. The projection of the matrix U = [u₁; …; u_s; …; u_S] can be expressed as a matrix V = [v₁; …; v_s; …; v_S], where the low dimensional hyperspectral vector v_s is the projection of u_s. The proposed algorithm is interpreted as Theorem 4.

Theorem 4.

For constants ε and β > 0, an integer K_PRP is chosen according to Equation (5). R_PRP is a D × K_PRP matrix whose entries belong to the standard normal Gaussian distribution. The N × K_PRP sub-matrix V^m is defined based on Equation (6) to present the projection of U^m in K_PRP dimensional subspaces. Then, for the same ε, β, K_PRP and R_PRP, with the probability of no distortion at least

P_PRP = min [P_PRP1, P_PRP2] = 1 − 2 (S/M)^−β

(11)

for any two hyperspectral vectors u_s and u_s′ of matrix U, the distance of the hyperspectral vectors v_s and v_s′ of matrix V remains as Equation (4).

The PRP algorithm can greatly reduce the dimensionality of HSIs, thereby reducing the computational load of subsequent HSI classification. The probability of no distortion P_PRP increases as the parameter β increases and the number of sub-matrices M decreases based on Equation (11). Theorem 4 demonstrates that the PRP algorithm does not introduce a significant distortion in HSIs. Meanwhile, the proposed algorithm can also gain a high probability of no distortion in the case of dividing a small number of sub-matrices. In this case, the proposed algorithm can approximate the distance between each hyperspectral vector pair while projecting HSI into a space with a lower dimensionality than the RP algorithm. Thus, considering Theorem 1 and Theorem 4 comprehensively, the PRP algorithm is feasible in calculation and dimensionality reduction for large-size HSI. Moreover, when the number of sub-matrices is too small, the proposed algorithm cannot ensure the achievement of the purpose of dimensionality reduction. Therefore, this paper gives the minimum number of sub-matrices that guarantees high dimensionality can be reduced. Additionally, it is defined as Lemma 1. (Proof of Lemma 1 can be found in Appendix B).

Lemma 1.

Suppose that 1.5 > ε > 0 and β > 0, then

$$M>S\times e^{-\frac{D}{24}}\iff K_0\le D$$

(12)

Figure 2 presents an example of dividing an HSI with five million hyperspectral vectors into one million sub-HSIs. The number of hyperspectral vectors N in each sub-HSI is five in Figure 2, which is calculated by dividing the number of hyperspectral vectors S in the HSI by the number of sub-HSIs M, that is, N = S/M = 5,000,000/1,000,000 = 5.

3.3. HSI Classification Based on the PRP Algorithm

The classification based on the PRP algorithm is introduced from two aspects. One is the optimization strategy of the projection matrix to increase the class separability of the dimensionality reduction result with the assistance of samples, and then the optimization strategy is applied to the PRP algorithm. The other is the HSI classification algorithm integrates the PRP algorithm for dimensionality reduction and MD classifier for classification.

3.3.1. Optimization Strategy of the Projection Matrix

Because the projection matrix is generated randomly, there is a case where the dimensionality reduction result with low class separability is not suitable for subsequent HSI classification. Therefore, it is necessary to exploit a projection matrix that is applicable for the classification task. The optimization strategy of the projection matrix is to increase the class separability of the low dimensional HSI. The class dissimilarity measurement is defined as large inter-class distance and small intra-class variance. The calculation method is the sum of the distance between classes divided by the variance within the class. Specifically, with the aid of samples, the projection matrix that maximizes the class dissimilarity is selected among multiple samplings.

Samples matrix for all classes of an HSI is expressed as X. Specifically, X = [X¹; …; X^l; …; X^L], where l is the index of classes and L is the number of classes known as a prior. And X^l is the lth samples matrix,

$${\mathit{X}}^l=\left[\begin{array}{cccc}{x}_{11}^l& x_{12}^l& \cdots & x_{1D}^l\\ x_{21}^l& x_{22}^l& \cdots & x_{2D}^l\\ ⋮& ⋮& \ddots & ⋮\\ x_{H1}^l& x_{H2}^l& \cdots & x_{HD}^l\end{array}\right]$$

(13)

where H is the number of samples in a single class, which is the same for each class. The selection of the projection matrix R_PRP is very important for the HSI classification. To be specific, many matrices are generated based on the standard normal Gaussian distribution. These matrices constitute a sampling matrix of the projection matrix, Q = [Q¹, …, Q^t, …, Q^T], where t is the index of sampling, T is the number of samplings, and Q^t is the projection matrix of the tth sampling. Then, each projection matrix Q^t is used to project the samples matrix X. The low dimensional samples matrix W = [W¹, …, W^t, …, W^T] is the projection of samples matrix X by the sampling matrices Q. Moreover, its calculation is derived from Equation (6), that is, W^t = XQ^t/√K_PRP. Specifically, W^t = [W^t1, …, W^tl, …, W^tL], where W^tl is the low dimensional samples matrix of class l of the tth sampling. That is,

$${\mathit{W}}^{tl}=\left[\begin{array}{cccc}{w}_{11}^{tl}& w_{12}^{tl}& \cdots & w_{1K_{\mathrm{PRP}}}^{tl}\\ w_{21}^{tl}& w_{22}^{tl}& \cdots & w_{2K_{\mathrm{PRP}}}^{tl}\\ ⋮& ⋮& \ddots & ⋮\\ w_{H1}^{tl}& w_{H2}^{tl}& \cdots & w_{H{K}_{\mathrm{PRP}}}^{tl}\end{array}\right]$$

(14)

The feature mean hyperspectral vector of low dimensional samples matrix W^tl of class l of the tth sampling is described as

$${\mathit{W}}_{\mathrm{mean}}^{tl}=\frac{1}{H}{\mathit{W}}^{tl}=\frac{1}{H}\left[\begin{array}{cccc}{w}_{11}^{tl}+\cdots +w_{H1}^{tl}& w_{12}^{tl}+\cdots +w_{H2}^{tl}& \cdots & w_{1K_{\mathrm{PRP}}}^{tl}+\cdots +w_{H{K}_{\mathrm{PRP}}}^{tl}\end{array}\right]$$

(15)

Next, distance between lth class and l′th class divided by the variance of the lth class is defined as follows,

$$I_{l{l}^{\prime }}^t=\frac{{‖{\mathit{W}}_{\mathrm{mean}}^{tl}-{\mathit{W}}_{\mathrm{mean}}^{t{l}^{\prime }}‖}^2}{var\left({\mathit{W}}_{\mathrm{mean}}^{tl}\right)}$$

(16)

where l′ is the index of classes. Then, the class dissimilarity is calculated as follows,

$$J^t={\displaystyle \sum _{l=1}^L{\displaystyle \sum _{l^{\prime }=1}^L{I}_{l{l}^{\prime }}^t}}$$

(17)

The better projection matrix is the t*th sampled projection matrix with the large class dissimilarity. Then,

$$t^{*}=\arg \underset{t=1,\dots ,T}{\max }{J}^t$$

(18)

Finally, the projection matrix R_PRP takes to the t*th sampled projection matrix Q^t*, where t* is obtained by Equation (18),

$${\mathit{R}}_{\mathrm{PRP}}={\mathit{Q}}^{t*}$$

(19)

The projection matrix selection strategy can increase the separability of dimensionality reduction results, thereby increasing the accuracy of HSI classification.

3.3.2. HSI Classification Based on Optimization Strategy of the Projection Matrix

HSI classification algorithm is described in this section. First, the low dimensional samples matrix Y = [Y¹; …; Y^l; …; Y^L] is the projection of samples matrix X by the projection matrix obtained according to Equation (19) with K_PRP dimensionalities. To be special,

$${\mathit{Y}}^l=\left[\begin{array}{cccc}{y}_{11}^l& y_{12}^l& \cdots & y_{1K_{\mathrm{PRP}}}^l\\ y_{21}^l& y_{22}^l& \cdots & y_{2K_{\mathrm{PRP}}}^l\\ ⋮& ⋮& \ddots & ⋮\\ y_{H1}^l& y_{H2}^l& \cdots & y_{H{K}_{\mathrm{PRP}}}^l\end{array}\right]$$

(20)

is the lth low dimensional samples matrix. The feature mean hyperspectral vector of the lth low dimensional samples is expressed as

$${\mathit{Y}}_{\mathrm{mean}}^l=\frac{1}{H}{\mathit{Y}}^l=\frac{1}{H}\left[\begin{array}{cccc}{y}_{11}^l+\cdots +y_{H1}^l& y_{12}^l+\cdots +y_{H2}^l& \cdots & y_{1K_{\mathrm{PRP}}}^l+\cdots +y_{H{K}_{\mathrm{PRP}}}^l\end{array}\right]$$

(21)

Consequently, an HSI is denoted by U = [u₁, …, u_s, …, u_S], where s is the index of hyperspectral vectors, S is the number of hyperspectral vectors, and u_s is a hyperspectral vector with D dimensionalities. The low dimensional HSI V = [v₁, …, v_s, …, v_S] is obtained by the PRP algorithm, where v_s is a low dimensional hyperspectral vector with K_PRP dimensionalities. Then, the final classes to which the low dimensional HSIs belong are determined by executing the MD classifier. That is to calculate the distance between the low dimensional hyperspectral vector and the samples feature mean hyperspectral vector of each class.

The distance matrix is characterized as Z = [z₁, …, z_s, …, z_S], where z_s is the sth distance vector. Especially, z_s = [z_s1, …, z_sl, …, z_sL], where z_sl is the distance between a low dimensional hyperspectral vector v_s and feature mean hyperspectral vector Y^l_mean according to Equation (21), and the calculation method is as follows,

z_sl = ||Y^l_mean − v_s||

(22)

The determined classification result can be defined as f = [f₁, …, f_s, …, f_S] in the HSI classification process. That is, the class with the smallest distance is the class to which the low dimensional hyperspectral vector belongs,

$$f_s=\arg \min {\mathit{z}}_s=\arg \underset{l=1,\dots ,L}{\min }\left\{z_{sl}\right\}$$

(23)

3.4. The Complexity of the Proposed Algorithm

The space and time complexity of the proposed algorithm are analyzed here. This section studies the complexity in two parts: the optimization strategy of the projection matrix and the HSI classification algorithm based on the optimization strategy of the projection matrix.

The main contribution of the complexity of the optimization strategy of the projection matrix is to calculate the projection matrix Q. To update Q, it takes O(HLD) space and O(HLDK_PRP) time for the calculation of the low dimensional sample matrix W. Furthermore, O(HK_PRP) space and O(THL²K_PRP) time are required to calculate the class dissimilarity J^t. The main contribution of the complexity of the HSI classification algorithm based on the optimization strategy of the projection matrix is to calculate the distance matrix Z. To update Z, it takes O(MND) space and O(MNDK_PRP) time for the calculation of the low dimensional image V. Furthermore, O(MNK_PRP) space and O(MNLK_PRP) time are required to calculate the distance z_sl. To sum up, the overall space complexity of the proposed algorithm is O(MND), and the overall time complexity of the proposed algorithm is O(HLDK_PRP + THL²K_PRP + MNDK_PRP + MNLK_PRP).

3.5. The Flow Chart of the Proposed Algorithm

The detailed process of the proposed Algorithm 1 is described as follows.

Algorithm 1 The detailed process of the proposed classification algorithm

Input: HSI U and samples X. Output: Classification result f. Step 1. Initializing parameters ε, β, L, T, H, M, S. Step 2. Dividing an HSI evenly U = [U1; …; UM]. Step 3. Calculating the lowest projection dimensionality K0PRP by using Equation (5). Step 4. Forming the projection matrix RPRP by using Equation (19). Step 5. Generating a low dimensional sub-HSI Vm by using Equation (6). Step 6. Getting the low dimensional HSI V = [V1; …; VM]. Step 7. Generating the feature mean hyperspectral vector of the low dimensional samples of lth class Ylmean by using Equation (21). Step 8. Calculating the distance zsl by using Equation (22).

In addition, Figure 3 illustrates a flow chart of the proposed classification algorithm in this paper.

4. Experiments and Results

To validate the effectiveness of the proposed algorithm, classification experiments in real HSIs were performed on a PC with Intel (R) Xeon (R) CPU E7-8880 v3, 2.30 GHz and 96.0 GB memory using MATLAB R2020a. The experimental parameters are set as shown in

Table 2. Projection parameter settings in the experiment.

	T	H	M	KPRP	KTRP	KRP
Pavia Centre image	10	10	36,598	33	100	349
Salinas image	10	10	2295	66	86	299
Chikusei image	10	10	3145	33	79	275
LongKou image	10	10	102,271	21	106	367

. The table contains the number of samplings T, the number of samples in each class H, the number of sub-HSI M, and low dimensionalities. It can be seen from Table 2 that the K_PRP is less than three-quarters of that of the K_TRP and one-eighth of that of the K_RP, which greatly reduces the amount of computation. In addition, the K_RP is larger than the high dimensionality D, which means that the RP algorithm cannot achieve the purpose of dimensionality reduction.

To verify the effectiveness of the projection matrix selection strategy, taking the LongKou image as an example, the randomly generated projection matrix (projection matrix without selection) is used as a comparison algorithm. Figure 4 shows the change of the spectral curve of the mean vector of samples of 9 classes of the LongKou image to evaluate the distance change between classes in the dimensionality reduction results. The horizontal coordinates represent the number of bands, and the vertical coordinates represent the intensity value. Figure 4a–c represent spectral curves of all classes without dimensionality reduction, dimensionality reduction without matrix selection, and dimensionality reduction based on matrix selection, with legends for all classes, respectively. Comparing the three images in Figure 4, the low dimensional spectral curves obtained by the projection matrix selection strategy have good separability for each class in each dimensionality, which fully reflect the superiority of the projection matrix optimization strategy.

For the sake of validating the proposed classification algorithm, the TRP-MIV algorithm, the algorithm for combining PRP and MGSR (PRP-MGSR), and the CAFCM algorithm are used as comparison algorithms. It is worth noting that the classification experiments of all experimental images have been tested 100 times. For all algorithms, the classification results with the highest accuracy in 100 trials are exhibited in Figure 5. Figure 5(a1–a4,b1–b4,c1–c4,d1–d4) show the classification results for four experimental images by the proposed algorithm and the TRP-MIV, PRP-MGSR, and CAFCM algorithms, respectively. Additionally, Figure 5(a5–d5) are legends for all classes of the four experimental images. Meanwhile, to better show the classification effect of the proposed algorithm, Figure 6 superimposes the final classification results of HSIs on the respective false-color images (see Figure 1(a1–d1) for details). Figure 6 shows the outlines of the superposition results of the four experimental images of the four algorithms.

From a qualitative point of view, compared with the other three comparison algorithms, the classification performance of the proposed algorithm is superior, which proves the reliability of the proposed classification algorithm for large-size HSIs. According to qualitative experimental results, the TRP-MIV algorithm has misclassification for complex HSIs with large hyperspectral vectors. In addition, comparing the PRP-MGSR algorithm, the proposed algorithm simultaneously considers the inter-class distance and intra-class variance and can achieve better classification results. The proposed algorithm preserves the class separability so that each class is classified exactly. In theory, the more dimensionalities are reduced, the more the loss of spectral information will be, which results in an unsatisfactory classification effect. However, the proposed algorithm attains a classification effect better than the CAFCM algorithm for four experimental images. Moreover, the CAFCM algorithm needs to be projected into a space with a higher dimensionality than the high dimensionality so that it cannot achieve the effect of dimensionality reduction.

The kappa coefficient, overall accuracy (OA), average accuracy (AA), average precision rate (APR), and running time are calculated to effectively evaluate the classification performance.

Table 3. Accuracy evaluation of classification results for the Pavia Centre image. (The values in brackets are the variance of 100 trials, and those not in are the mean).

	The Proposed Algorithm	TRP-MIV	PRP-MGSR	CAFCM
Kappa coefficient	0.83 (0.02)	0.81 (0.05)	0.79 (0.02)	0.38 (0.11)
OA/%	89.65 (1.03)	88.02 (2.94)	87.34 (1.25)	53.72 (11.68)
AA/%	78.99 (1.94)	75.43 (5.88)	73.22 (3.02)	38.03 (5.65)
APR/%	79.80 (1.70)	74.49 (5.81)	72.72 (2.60)	34.11 (6.01)
Running time/s	1.16 (0.09)	2.93 (0.29)	10,209.72 (1538.24)	5798.18 (855.03)

,

Table 4. Accuracy evaluation of classification results for the Salinas image. (The values in brackets are the variance of 100 trials, and those not in are the mean).

	The Proposed Algorithm	TRP-MIV	PRP-MGSR	CAFCM
Kappa coefficient	0.91 (0.00)	0.87 (0.01)	0.88 (0.00)	0.59 (0.14)
OA/%	91.98 (0.23)	89.00 (0.71)	90.20 (0.38)	63.77 (11.47)
AA/%	85.97 (0.92)	84.24 (1.25)	80.14 (1.70)	48.95 (11.68)
APR/%	79.89 (0.66)	75.57 (1.20)	75.28 (1.33)	45.32 (17.93)
Running time/s	0.37 (0.10)	3.93 (0.11)	100.08 (8.65)	1623.16 (226.78)

,

Table 5. Accuracy evaluation of classification results for the Chikusei image. (The values in brackets are the variance of 100 trials, and those not in are the mean).

	The Proposed Algorithm	TRP-MIV	PRP-MGSR	CAFCM
Kappa coefficient	0.97 (0.01)	0.92 (0.05)	0.90 (0.01)	0.35 (0.10)
OA/%	97.74 (0.38)	94.15 (3.36)	92.93 (0.74)	48.52 (7.77)
AA/%	96.66 (1.06)	87.89 (8.85)	91.63 (0.81)	36.58 (9.75)
APR/%	94.57 (0.56)	87.29 (7.67)	85.00 (1.93)	33.62 (6.54)
Running time/s	0.21 (0.76)	1.59 (0.16)	52.20 (11.34)	289.21 (34.70)

and

Table 6. Accuracy evaluation of classification results for the LongKou image. (The values in brackets are the variance of 100 trials, and those not in are the mean).

	The Proposed Algorithm	TRP-MIV	PRP-MGSR	CAFCM
Kappa coefficient	0.82 (0.01)	0.62 (0.01)	0.62 (0.03)	0.47 (0.01)
OA/%	86.24 (0.51)	68.90 (0.81)	69.73 (2.93)	56.84 (1.03)
AA/%	77.94 (1.57)	56.59 (2.12)	55.51 (3.51)	33.69 (2.56)
APR/%	74.23 (1.50)	52.05 (1.65)	52.78 (2.36)	39.04 (3.21)
Running time/s	2.42 (0.24)	7.76 (0.40)	7215.75 (835.60)	16,698.02 (1512.06)

reveal the quantitative evaluation results of four experimental images, respectively. The values in the tables are the mean and variance of the kappa coefficient, OA, AA, APR, and running time calculated by different algorithms for 100 trials, in which the values in brackets are the variance and those not in are the mean.

From a quantitative perspective, the classification results of the proposed algorithm surpass other algorithms. The running time of the proposed method for four image classification experiments is about 1.16, 0.37, 0.21, and 2.42 s, and the mean values of their kappa coefficient are larger than 0.82. All evaluation parameters intuitively state that the proposed algorithm accurately classifies the real HSIs. According to quantitative experimental results, although the mean values of AA and APR are less than 80% for the Pavia Centre image, the mean value of OA is more than 89%. For the Salinas image, the variance value of the kappa coefficient of the proposed algorithm in 100 trials is approximately equal to zero, which indicates that the classification results of the proposed algorithm are relatively stable. In addition, for the Chikusei image, the mean values of all evaluation parameters of the proposed algorithm are at least 0.05, 3.59%, 4.73%, 7.28%, and 7.57 times larger than those of the comparison algorithms, which fully reflects the superiority of the proposed classification algorithm. Although the kappa coefficient of the proposed algorithm is not very high for the LongKou image, the accuracy of the proposed algorithm is still higher than that of all comparison algorithms. Additionally, the values of OA, AA, and APR outperform the comparison algorithms by at least 16.51%, 21.35%, and 35.19%. In summary, the proposed algorithm can obtain reliable classification results in a very short time.

5. Discussion

In order to verify the superiority of the proposed classification algorithm, all classification experiments have been tested 100 times. The TRP-MIV algorithm [38] has better classification results for HSIs with a large inter-class distance. The PRP-MGSR algorithm is relatively stable in classification results in multiple tests. This is because the MGSR algorithm [14] has strong applicability to various data, but the running time is relatively long. The CAFCM algorithm [39] is more likely to be disturbed by noise, and the classification results of the algorithm are better for images with less noise. The proposed classification algorithm considers the class dissimilarity, which can reach the ideal classification result in a shorter time.

In order to use the proposed algorithm more widely, this section discusses the influence of the projection parameters on dimensionality reduction and the influence of the number of samples and the number of samplings on the classification accuracy. First, the focus is to analyze the effects of parameters ε, β, and the number of sub-HSIs M on the lowest projection dimensionality K⁰_PRP. According to Equation (5), taking the Salinas image as an example, Figure 7 represents the change curves of the lowest projection dimensionality with different parameters, where the horizontal coordinates are different parameters, and the vertical coordinates are K⁰_PRP. Figure 7a is the change curve of ε and K⁰_PRP when β = 0.5 and M = 2295 (according to Table 2), where the position of the black dotted line is the inflection point. When ε is equal to 1.5, the denominator term of Equation (5) is 0, so the curve in Figure 7a is discontinuous at ε = 1.5. When ε < 1.5, K⁰_PRP first decreases and then increases with the increase of ε and takes the minimum value when ε = 1. When ε > 1.5, K⁰_PRP increases with the increase of ε and is always negative, which is meaningless for dimensionality reduction. Figure 7b is the change curve of β and K⁰_PRP when ε = 1 and M = 2295. As seen in Figure 7b, K⁰_PRP increases with the increase of β. Figure 7c is the change curve of M and K⁰_PRP when ε = 1 and β = 0.5. The position of the black dotted line is the number of bands, and the values on the upper side and the left side of the boundary line cannot achieve dimensionality reduction. As seen from Figure 7c, K⁰_PRP decreases with the increase of M. Therefore, to reduce computation and memory, it is necessary to project the HSI to a smaller dimensional space. Then, it is best to set ε to 1 within (0, 1.5); β should be set as small as possible, and M should be set as large as possible. Next, in order to reflect the impact of the number of samples H and the number of samplings T on classification performance, the change curves in the values of OA with different T and H are shown in Figure 8. The number of samples and samplings are taken from 10 to 650, respectively. Figure 8 shows that as T and H increase, the OA value varies only in a small range. The larger the number of samplings and samples, the smaller the range of variation. Among them, Figure 8a shows that with the increase of T, the error of the OA value is 1%, and Figure 8b shows that with the increase of H, the error of the OA value is 2.5%. Therefore, in practical applications, only a small number of samples and few projection matrix samplings can be used to obtain better classification results. In addition, in order to ensure that the PRP algorithm can still maintain the distance after the HSI is divided, the number of samples in each class in this paper is set to be the same. The focus of future research will be on dividing an HSI into sub-HSIs with different numbers of hyperspectral vectors to improve utility.

6. Conclusions

A new dimensionality reduction method, the PRP algorithm, is presented for large-size HSIs in this paper. It is worth mentioning that this paper also theoretically proves the distance preservation property of the PRP algorithm in detail. The PRP algorithm projects the large-size HSI into a space with a lower dimensionality than the RP algorithm. Therefore, when applied to the large-size HSI, the PRP algorithm no longer has the problem of the inability to reduce the dimensionality. In addition, the proposed algorithm utilizes the large inter-class distance and small intra-class variance as the class dissimilarity measurement for selecting the projection matrix, which can preserve the class separability of the low dimensional HSI well. Compared with TRP-MIV, PRP-MGSR, and CAFCM algorithms, the proposed algorithm has the shortest running time and the highest classification accuracy.