Feature Extraction of Hyperspectral Images Based on Subspace Band Selection and Transform-Domain Recursive Filtering

2.1 Lasso-based band selection

For hyperspectral image data, there is a closer correlation between adjacent bands than non-adjacent bands. Thus, it is entirely feasible to sort the bands of a hyperspectral image using the Lasso algorithm. The steps of the algorithm are as follows:

For a hyperspectral image containing l types of land covers and N spectral bands, M training samples are given. Each sample has a known spectral type and sample label. Then, each pair of samples can be expressed as (x_i, y_i), where i \in[1,2, \cdots, N]. Sample x_i is an M-dimensional column vector of the spectral response of the i-th sample. Obviously, the j-th element of x_i is the spectral response of the j-th band of the hyperspectral image; y_i is an l-dimensional column vector, representing the class labels of all i samples.

The Lasso algorithm introduces the l₁ norm as a penalty term to the ordinary least squares (OLS) estimation, and ensures the sparsity of the solution of the linear regression problem. Then, the sparsity coefficient of each band can be obtained by:

$\hat{\alpha}=\underset{\alpha}{\arg \min } \sum_{i=1}^{N} \frac{1}{2}\left\|\boldsymbol{y}_{i}-\alpha^{T} \boldsymbol{x}_{i}\right\|_{2}^{2}+\lambda\|\alpha\|_{1}$      (1)

where, y_i is the class label of a sample; coefficient α is an N+1-dimensional column vector; N is the dimensions of the original hyperspectral image (i.e., the number of bands).

Since the norm is singular on the vertex of each orthonormal basis, the coefficient α calculated by formula (1) is sparse. λ is the regularization parameter used to adjust the sparsity of the solution.

2.2 Transform-domain recursive filtering

Transform-domain recursive filtering is an edge-preserving filter. It is capable of preserving the edges, lines and other details of the image, while smoothing the texture and noise in the image. Thanks to this advantage, the transform-domain filter and other edge-preserving filters have piqued much interest among researchers in the field of image processing. For example, Wang et al. [16] used edge-preserving filters to extract features from hyperspectral images, and fully utilized the spectral information and spatial information of the image in the post-processing of pixel-based classification results. Luo et al. [17] proposed a spectral-spatial feature extraction method based on edge-preserving filtering. The method removes unimportant parts of hyperspectral images, and further improves the interpretability of hyperspectral image classification results. Zhu et al. [18] decomposes the hyperspectral image through independent component analysis (ICA), and performs edge-preserving filtering on image edges, in order to extract the features for image classification.

This section only briefly describes the transform-domain recursive filter. The complete description is available in Gastal’s work [19].

Suppose there is a one-dimensional (1D) input signal I. Following the approximate distance-preserving transformation method, the input signal can be transformed to the transform domain Φ:

$\boldsymbol{V}_{i}=\boldsymbol{I}_{0}+\sum_{k=1}^{i}\left(1+\frac{\sigma_{s}}{\sigma_{r}}\left|\boldsymbol{I}_{j}-\boldsymbol{I}_{j-1}\right|\right)$     (2)

where, V_i, originally a function of distance calculation, is the transform domain signal; σ_s is the standard deviation of the spatial domain, responsible for adjusting the window size of the filter; σ_r is the standard deviation of the value domain, responsible for adjusting the fuzziness of the filter. Then, the input signal I is subject to recursive filtering:

$\boldsymbol{W}_{i}=\left(1-\omega^{\rho}\right) \boldsymbol{I}_{i}+\omega^{\rho} \boldsymbol{W}_{i-1}$     (3)

where, W_i is the filter output for the i-th pixel in the input signal I; $\rho=\exp \left(-\frac{\sqrt{2}}{\sigma_{s}}\right)$, $\rho \in[0,1]$ is a feedback signal; ω is a distance coefficient representing the distance between two adjacent pixels W_i and W_i-1 in the transform domain. Obviously, as the feedback signal increases, ω^ρ will gradually decrease until reaching zero. Then, the propagation chain will be stopped, and the edges of the signal will be preserved.

For two-dimensional (2D) images, transform-domain recursive filter processes the image by performing 1D operations along each dimension of the image. Sun and Du [20] demonstrated that the iterative execution of three 1D filtering operations can obtain filtered images without artifacts.

4.1 Experimental setup

(1) Datasets

To verify the classification performance of BSTDRF on hyperspectral images, two real hyperspectral datasets were adopted for experiments: the Indian Pines dataset (Figure 2) and the University of Pavia dataset (Figure 3).

2.png

Figure 2. Indian Pine dataset (a) Pseudo-color image of the Indian Pines image, (b) Reference data, (c) Class names

3.png

Figure 3. University of Pavia dataset (a) Pseudo-color image of the University of Pavia image, (b) Reference data, (c) Class names

In 1992, the Indian Pines image was taken by an airborne visible infra-red imaging spectrometer (AVIRIS) in the Indian Pines lab base in northwestern Indiana, United States. This remote sensing image (size: 145 × 145 × 22) has a spatial resolution of 20m, and a wavelength range of 0.4-2.5μm. Before the hyperspectral image classification experiment, the spectral bands covering water absorbing areas must be removed from the original hyperspectral data: [104-108], [150-163], and 220. The removal reduced the number of spectral bands from 220 to 200.

The Indian Pines image contains multiple native vegetation, including farmland, forest, and others. According to the current exploration, the scene can be divided into 16 different targets: Class 1 represents Alfalfa, Class 2 represents Corn-no till, Class 3 represents Corn-min till, Class 4 represents Corn, Class 5 represents Grass/pasture, Class 6 represents Grass/trees, Class 7 represents Grass/pasture-moved, Class 8 represents Hay-windrowed, Class 9 represents Oats, Class 10 represents Soybeans-no till, Class 11 represents Soybeans-min till, Class 12 represents Soybeans-clean till, Class 13 represents Wheat, Class 14 represents Woods, Class 14 represents Buildings, and Class 16 represents Stone-steel towers.

The University of Pavia image was collected by a German reflective optics system imaging spectrometer from a city near the University of Pavia, Italy. This airborne image has a spatial resolution of 1.3m, and a wavelength range of 0.43-0.86μm. The 610 × 340 sized image contains 115 bands. Before the experiment, the 12 most noisy bands were removed.

There are nine different ground covers in the scene of University of Pavia: Asphalt (Class 1), Meadows (Class 2), Gravel (Class 3), Trees (Class 4), Metal sheets (Class 5), Soil (Class 6), Bitumen (Class 7), Bricks (Class 8), and Shadows (Class 9).

(2) Quality metrics

After classifying each hyperspectral image, it is necessary to evaluate the classification effect. The classification accuracy should be assessed against ground reference data. This paper adopts three common metrics of the classification accuracy of hyperspectral images, namely, overall accuracy (OA), average accuracy (AA), and Kappa coefficient. Among them, OA measures the percentage of pixels being classified correctly; AA refers to the mean percentage of correctly classified pixels in each class; Kappa coefficient excels in estimating the recognition accuracy, in the light of the impact of uncertainties on the classification results.

4.2 Classification results

(1) Parameter analysis

4-1.png

(a) Influence of σ_s

4-2.png

(b) Influence of σ_r

Figure 4. Analysis of the parameters σ_r and σ_s on Indian Pines data set

The proposed BSTDRF adopts the transform-domain recursive filtering technique. It is necessary to understand the roles of filter parameters σ_s and σ_r in classification, and determine their optimal values. That is why the Indian Pines image was selected for experiment.

During the experiment, 10% of samples from the Indian Pines dataset were chosen randomly as the training set, and the remaining 90% as the test set. Through cross validation, an empirical value was configured for σ_r. Different σ_s values were selected until the output was nearly optimal. Then, the σ_s value at the moment was recorded, and taken as the optimal filter parameter of BSTDRF. Next, the σ_s was fixed at the optimal value obtained in the first step, and different σ_r were selected until the output was nearly optimal. In this way, the two parameters of the transform-domain recursive filter were determined. The specific process is as follows:

Firstly, σ_r was fixed at 0.4, while σ_s was gradually increased from 10 to 100. The classification results are displayed in Figure 4(a). It can be learned that OA and AA reached the optimal level, when σ_s=70. After σ_s surpassed 70, OA, AA and Kappa coefficient all started to decline. This means an excessively large window of the filter will cause the loss of details. The spatial coefficient should be controlled within the effective range.

Secondly, σ_s was fixed at 70, while σ_r was gradually increased from 0.1 to 1.0. The classification results are displayed in Figure 4(b). It can be observed that all three indices fell on a low level, when σ_r was small. The undesired classification effect is attributable to the fact that a small filter parameter can only extract a limited number of features from a small neighborhood. When σ_r was large, all three metrics started to decrease. Thus, when the value domain parameter is large, the transform-domain recursive filter is approximately equivalent to a Gaussian filter, which makes the image fuzzy. When σ_r=0.4, the overall performance of the three parameters was the best. Overall, this paper sets σ_s as 70, and σ_r as 0.4.

(2) Comparison of different classification methods

Our BSTDRF was compared with seven widely used classifiers of hyperspectral images, namely, SVM, PCA-based SVM, ICA-based SVM, multidimensional scaling (MDS)-based SVM [21], extended morphological profiles (EMP) [22], linear regression for machine learning (LRML) [23], and DTB [24].

The SVM methods were realized using the LIBSVM library and the radial basis function (RBF) kernel. The parameters of each SVM classifier were determined through five-fold cross-validation. All experiments were conducted using MATLAB 2019 on a computer with Intel Core i7-8700 CPU and 16 GB RAM. To ensure the objectivity of the evaluation, each experiment was repeated 10 times with random training samples. For PCA and ICA, the top-20 principal or independent components were imported to the filter, and the default parameter setting was adopted for the remaining methods.

The first experiment was carried out on the Indian Pines dataset. 10% of the samples were randomly selected as training samples, and the remaining 90% were used as test samples. Figure 5 shows the classification results and OAs of different algorithms. It can be observed that the EMP achieved better classification accuracy (OA=89.1%) than pixel-level feature extraction algorithms, by utilizing the spatial structure information in the image. However, some noises remained in the classification results of EMP. LRML, DTB and BSTDRF all effectively removed the incorrect classification, which looks like noises, and achieved better recognition accuracy. Our BSTDRF achieved the best OA, which is 20.2% higher than that of SVM, 5.1% higher than that of LRML, and 0.4% higher than that of DTB. Hence, reconstructing the hyperspectral image with important bands can reduce the spectral dimensionality of the hyperspectral image, eliminate noisy bands, and retrain the useful information in the images.

5.png

Figure 5. Classification results obtained by different methods on Indian Pines dataset (a) SVM, (b) PCA-SVM, (c) ICA-SVM, (d) MDS-SVM, (e) EMP, (f) LRML, (g) DTB, (h) Proposed BSTDRF

Table 1 displays the OA, AA and Kappa coefficient of each method on each class in the scene. It can be seen that, when the training sample took up 10% of ground reference data, the SVM could effectively differentiate between ground targets with very different spectrums, such as Grass_M, Grass_T, Grass_P, Hay_W, Wheat, Woods, and Stone. However, this pixelwise classifier failed to recognize ground targets with similar spectrums, due to the neglection of the spatial information of the hyperspectral image. For instance, the classification accuracy on Corn was merely 43.7%, and that on Buildings was only 53.9%.

Compared with the classification directly on the original data of the hyperspectral image, PCA and ICA could effectively reduce the dimensionality of the data, but their OA declined by 2.3% - 4.6%.

EMP relies on the multiscale morphological operator to extract the spatial features from the hyperspectral image for classification. LRML introduces multilayer regression into the spatial information of the image, and optimizes the spectral classification result of logistics regression classifier. DTB introduces the spatial information of the image via multi-filter fusion, and thus optimizes the space of the spectral classification results. As shown in Table 1, these three spatial-spectral feature extractors outperformed the pixelwise classifier.

The proposed BSTDRF realized the highest OA, AA and Kappa coefficient. Compared with the original SVM, our approach improved the OA from 79.2% to 95.2%, the AA from 78.3% to 94.6%, and the Kappa coefficient from 76.9% to 95.4%. Some ground targets, namely, Buildings and Alfalfa, cannot be identified correctly by some pixelwise classifiers. Our approach can improve their classification accuracy to 91.5% and 81.3%, respectively.

The experiment on the Indian Pines scene suggests that our BSTDRF has better classification performance than other classification methods.

Table 1. Classification accuracies of different methods on the Indian Pines scene (10% for training set)

6.png

Figure 6. Classification results obtained by different methods on University of Pavia data set (a) SVM, (b) PCA-SVM, (c) ICA-SVM, (d) MDS-SVM, (e) EMP, (f) LRML, (g) DTB, (h) Proposed BSTDRF

Table 2. Classification accuracies of different methods on the University of Pavia scene (4% for training set)

The second experiment was carried out on the University of Pavia dataset, which has many more ground references than Indian Pines dataset. Hence, 4% of samples were selected randomly for the training set, and the rest 96% for the test set. Figure 6 shows the classification results and OAs of different algorithms.

It can be seen that the four spatial-spectral classifiers, including EMP, LRML, DTB, and our approach, achieved better classification accuracy than the pixel-wise classification methods. In addition, the OA (96.2%) of our approach was better than that of the other three spatial-spectral classifiers. The advantage was 1.8% over EMP, 1.5% over LRML, and 1.2% over DTB. This is because our approach adopts the Lasso-based band selection strategy, which reduces the dimensionality of data, while preserving the physical meaning of the data.

Table 2 displays the OA, AA and Kappa coefficient of each method on each class in the University of Pavia scene. It can be seen that, among the four pixel-wise classifiers, the approach that directly imports the original data of hyperspectral images to the SVM, boasted the best effect. The other three pixel-wise methods, namely, PCA, ICA, and MDS, lost some useful information during dimensionality reduction, which affects the classification accuracy, although the dimensionality reduction improves calculation efficiency.

Taking the Gravel ground targets for example, the OA obtained after the dimensionality reduction by PCA, ICA, and MDS was 2.1%, 5.3%, and 1.5% lower than the OA achieved using the original data directly, respectively.

In addition, the pixel-wise classifiers fail to consider the spatial information of images, and thus performed poorly in recognizing some ground targets. For example, the highest OA of the pixel-wise classifiers was 69.8% (MDS), while the OA of our approach was 37.4% higher (95.9%).

The second experiment shows that our approach removes the noisy bands from the hyperspectral image through band selection, and retains the contours and edges of the image via transform-domain recursive filtering. Thus, our approach is suitable for precision agriculture, as well as urban planning.