3.1.

Model of Spectral Window

Given a hyperspectral feature vector $x={(x_1,x_2,\cdots ,x_D)}^T$, each component $x_d$ denotes a spectral reflectance or radiance value measured at the band $d$, $d=1,2,\dots ,D$. $D$ is the total number of spectral bands in a spectrometer. This number could be as many as 200 for a normal hyperspectral sensor. Given that $T_{w,b}$ is a transform that can choose a subset of contiguous bands at the spectrum position $b$ within a $w$-width window, $T_{w,b}(x)$ is a chunk of spectral signal within the spectral window (see Fig. 2).

To automatically find the best spectral window, we may consider MI as a metric $M$. As long as the derivative of this MI is known, we can decide how the spectral window should move, i.e., the gradient ascending algorithm. For other similarity metrics, such as correlation or spectral angle, it is not so easy to find their analytic expression and apply them to the gradient-based algorithms. According to Eq. (3), the MI can be written as

We can see that the MI described in the right side of the first row of Eq. (10) has two terms. The first one is the entropy of the class label variable $r$. It is not a function of the window-shifting transform $T$. The second term is a conditional entropy of $r$ given the chunked spectral feature vector $x_{w,b}$. When the class label variable $r$ and the chunked spectral feature vector $x_{w,b}$ are related, the amount of entropy, $H(r|x_{w,b})$, will reduce. For example, if $r$ can be predicted by $x_{w,b}$ reliably, then $r$ will become less uncertain when we have the observation of $x_{w,b}$. If we intuitively understand the entropy as the amount of uncertainty, the decrease of uncertainty means less of entropy. In other words, if $x_{w,b}$ is a good observation or informative feature subset to predict the unknown variable $r$, the uncertainty of $r|x_{w,b}$ will reduce more than other less discriminatory feature subsets. As a result, $H(r|x_{w,b})$ will decrease. From Eq. (10), the less amount of $H(r|x_{w,b})$ will increase the MI, $I(r,x_{w,b})$. Consequently, considering that $x_{w,b}=T_{w,b}(x)$, maximizing the MI in Eq. (10) will encourage the spectral window to shift to a wavelength region in which the spectral features would have a better capability to predict the class label.

Also, there is a bound relationship between the MI and the Bayes error,²⁵ so the bands subset $x_{w,b_0}$ found by Eq. (8) would have a good discriminatory capability for classification in the sense of Bayes error. In the following paragraphs, we discuss how to obtain the differentiable MI to maximize $I(r,x_{w,b})$.

3.2.

Differentiable Mutual Information

In Eq. (10), the MI is defined as the difference of three entropies, $H(r)$, $H(x_{w,b})$, and $H(r,x_{w,b})$. From Eq. (2), these entropies are defined in terms of sums over the probability densities associated with the random variable $r$ and $x_{w,b}$. In hyperspectral image classification, these densities are usually unavailable and need to be estimated. Considering that the estimated signal could be quite a high-dimensional ($>200$) variable, it is difficult to collect enough samples to validate its histogram. However, it is always possible to manually label a small number of training samples, such as 50 to 100 pixels, for a specific task. Based on this small number of samples, Parzen window method²⁶ can be applied to approximate the underlying probability density, which is briefly presented as follows.

Let the sample set of a random variable be $Y=\{y_1,y_2,\cdots ,y_N\}$; then the probability density function of $Y$ is estimated as the sum of a group of normalized window functions centered on the samples.

After modeling the probability density functions, we can further calculate the entropy of random variables by using the approximation formula in Eq. (11). This can be written as follows:

Using Eqs. (11) to (13), the analytic expression of the density functions and entropy can be derived, and the derivative of the entropy with respect to the transform $T$ can be obtained as follows:^23,27

Combining Eqs. (13), (14), and (15), the derivative of the entropy is given by

Since a gradient ascending approach is used to find the maxima of MI, it is necessary to calculate the derivative of the MI with respect to the window-shifting transform $T$. From Eq. (10), the derivative of $I(r,x_{w,b})$ can be represented by

Substituting the derivatives in Eq. (17) with the results in Eq. (16) [note that the $y$ in Eq. (16) could be a vector], the derivative of $I(r,x_{w,b})$ is calculated by

The derivatives $(d/dT)(w_j-w_i)$ and $(d/dT)({\mathbf{v}}_j-{\mathbf{v}}_i)$ can be approximated by the transform difference $\mathrm{\Delta }(w_j-w_i)$ and $\mathrm{\Delta }({\mathbf{v}}_j-{\mathbf{v}}_i)$. Since in this research the transform $T$ is used to model a behavior of spectral windows moving along the spectral axis, we can calculate the difference by subtracting the values of the current spectral window from the values of the spectral window shifted by one spectral band, i.e., in Eqs. (19) and (20).

In this section, by following the methodology introduced in the literature,^23,27 we reformulated the high-dimensional MI for our hyperspectral application. This work is based on the two-dimensional scenario originally developed for medical image registration in Ref. 27. We also deduced a new difference transform for the applied spectral window’s moving model [see Eqs. (19) and (20)]. To apply the differentiable MI to our hyperspectral research, we still need a new search strategy, which will be discussed in the next section.

3.3.

Searching Algorithm

After approximating the derivative of MI, the maxima of MI are found by the gradient ascending algorithm, which is detailed in the following steps:

The local maxima are found by repeating the steps 2 and 3 until convergence is detected or a fixed number of iteration times is reached.^23,27

To avoid the local maxima, a random-restart strategy is adopted, which runs an outer loop over the above gradient ascending in several wavelength regions. Every region can be obtained by dividing the whole spectrum evenly. Each outer loop chooses a random initial position in the region to start gradient ascending. If a new run of gradient ascending produces a better result than the stored one, it replaces the stored solution. After running out each region, the best final solution is found as the global submaxima. This strategy is simple but effective, which is evidenced by the experimental results shown in Sec. 4.

4.1.

Searching Spectral Window to Build a Reliable Reference Map

In the experiment, we first examine if the proposed scheme can find the optimal subset of bands to approximate the ground truth. We implement the algorithm introduced in Sec. 3 to search the optimal spectral window, with window widths varying from 20 to 50 bands, respectively. The number of samples for density estimation, i.e., $N$, is 50, and the number of samples for entropy estimation, i.e., $L$, is also 50. Selection of these two numbers is empirical in this research, and they are chosen to fit the specific application. Apparently, to estimate different probably density functions (PDFs), different $N$ and $L$ may be needed, and their values (i.e., the size of samples) should match the complexities of the approximated PDFs. Since we expect that the subset bands contain the most discriminatory information, they should achieve the highest classification accuracy of any other subsets. Hence, we compare the classification accuracies based on the bands within and outside the found spectral windows respectively. In the experiments, we chose support vector machines (SVMs)^29,30 as the classifiers due to the higher dimensionality of input data. Also, previous works applying SVMs to hyperspectral data classification have shown competitive performance with the best available classification algorithms.^28,31 However, other classification algorithms, including the unsupervised approaches, are also applicable since the band selection proposed in this research is independent of the classifier applied (actually, it can be categorized as a filtering method in the feature selection family).

Figures 3(a) to 3(d) present the accuracy results under window widths of 20, 30, 40, and 50 bands, respectively. The solid lines denote the classification accuracy achieved by the found spectral window, and the star points are the classification results based on the other 200 random sampled spectral windows. The dashed lines stand for the means of the classification results from the 200 random sampled windows. The positions of the sampled spectral windows are decided by random numbers generated by the machine. The bands in these spectral windows are contiguous (see Fig. 5). The number of sampling is initially chosen as 200, which is in favor of the narrow window-width, such as 10 bands. In the case of wide window-width, such as 30, 40, or 50 bands, several repeats of sampling will happen, but this will not affect our comparison. It is seen that the accuracies based on the found spectral windows are indeed on the top of the results of all 200 random sampled windows, which justified the effectiveness of the searching algorithm. In other words, maximizing MI can indeed lead to finding of a good estimation of ground truth (i.e., the ideal classification result).

Fig. 3

Comparison of classification accuracy: - stands for the accuracy achieved by the found spectral windows, and ∗ the 200 random sampled windows as a function of the random sample sequence order; the widths of spectral window are (a) 20 bands, (b) 30 bands, (c) 40 bands, and (d) 50 bands.

Therefore in this experiment, the window-width is chosen as 35 bands. The above experiments are mainly used for justification of the algorithm. In practical applications, we can use MI value as the indication of the classification accuracy and do not need to implement the real classification test.

The window-width $w$ is an application-related parameter. How to choose this parameter is crucial to this method. Basically, we can set the window-width by two approaches. The first one is by a priori knowledge. For example, if the classes to be categorized are known, we can refer to the existing spectral libraries and measure the widths of the absorption valleys (i.e., a particular spectrum where the light is absorbed by the constituent atoms or molecules; usually these spectra are like valleys in a spectral curve). Then, the size (in wavelength) of the spectral window can be chosen as the maximum of the absorption valleys’ widths. The second way is by a validation test using the training data. Figure 4(a) shows the maximal MI found by different window-widths. It is seen that when $w$ is set between 30 and 40 bands, a higher MI value is found. Therefore in this experiment, the window-width is chosen as 35 bands. Figure 4(b) illustrates the accuracy change with respect to the positions of spectral windows. In this case, the starting band number of this spectral window is band no. 14. It should be noted that the above experiments are mainly used for justification of the algorithm. In practical applications, we can use MI value as an indication of the classification accuracy and do not need to implement the real classification test.

Fig. 4

Choosing the (a) width and (b) validating the position of the spectral window.

We also verified the proposed method by visual comparison. In Fig. 5, we illustrate the spectral window found by the proposed method, together with three other random sampled spectral windows (to help visual observation, all hyperspectral images shown in this paper are transformed using a suitable false color palette). To facilitate this comparison, the width of the spectral window is set up as 10 bands. Figures 5(a) to 5(c) illustrate $3\times 10$ spectral bands corresponding to three 10-bands width random sampled windows. Figure 5(d) displays 10 band images corresponding to the spectral window found by the searching algorithm. Comparing these spectral images with the accompanied reference map in Fig. 5(e), it can be seen that the images found by the proposed method contain relatively more discriminatory information than other random sampled bands. In Fig. 5(d), we can see that the inherent scale varies across the bands and can roughly distinguish the outline of the reference map, whereas in the bands of Figs. 5(a) to 5(c), this becomes less clear or impossible (for example, the third random sampled window that is unfortunately located in the atmospheric water absorption area). The visual comparison is consistent with the numerical analysis in Fig. 3.

Fig. 5

Visual comparison of three random sampled spectral window with the optimal spectral window; window-width is 10 bands; (a) random sample window 1, (b) random sample window 2, (c) random sample window 3, (d) optimal window, and (e) the reference map.

Through the above accuracy comparison and visual observation, it is verified that the proposed algorithm can find a spectral window to maximize the classification accuracy. Based on this spectral window, we are able to construct the estimated reference map by averaging the bands within the spectral window, such as in Eq. (6). If a spectral window works well for the classification, it gives a good indication to show that the characteristics of these bands (e.g., the clustering or separability of pixels’ values) have good capability to predict the class label. This can give information similar (but in different coding) to the one that the ground truth can provide. So a spectral window that works well for classification will also be good to build a reference map. Meanwhile, because the pixel values of a band can be considered as another kind of coding of ground truth, the bands within the spectral window can be regarded as a series of weaker reference maps. The proposed method makes these bands positively correlated (see the window constraint), with higher classification accuracy (see MI evaluation and the searching strategy); averaging them can provide a better estimation of reference map by reducing noise and uncertainty. This is the simplest approach to build the reference map by taking advantage of the proposed method. After constructing the reference map, band selection can be carried out, which is described as follows.

4.2.

Results on Band Selection

The main objective of band selection is to remove redundant spectral bands without degrading the classification accuracy too much. The experiment was designed to assess the change of classification accuracy as spectral bands are progressively removed according to the ranked MI values. To compare with the previous research,^28,32,33 we use a subscene comprising of four classes: corn-notill, soybeans-notill, soybeans-min, and grass/trees. Figure 6 shows the results for the three cases, where the MI was calculated with respect to the reference map accompanied with the dataset (the solid curve marked with triangle), the subjective estimated reference map (the dotted curve marked with circle) (note that this map is estimated by visual inspection from domain expert), and the optimally constructed reference map using the approach introduced in Sec. 3 (the dashed curve marked with rectangle), respectively. Data points in the figure are at 5% increments corresponding to removal of 11 bands at each step. The results depicted in Fig. 6 show that in this particular application scenario,

Fig. 6

Comparison of overall classification accuracy, AVIRIS 92AV3C scene.

Fig. 7

Individual mean classification accuracies on the four different classes: (a) corn-notill, (b) grass/trees, (c) soybeans-min, (d) soybeans-notill.

The results summarized in the above 2 showed another improvement was achieved by the proposed method: monotonic performance degradation. This is an expected result and means that the ranking order for bands-cutting is more consistent with their contribution to the classification accuracy. In contrast, the performance curves of the other benchmarked methods showed no such result. Finally, the accuracy degradation of the proposed method is slower than other methods, indicating a better robust performance for bands-cutting.

To justify the proposed method furthermore, we carried out an experiment based on a high-resolution hyperspectral dataset, i.e., Pavia University scene of ROSIS sensor (by courtesy of Prof. Paolo Gamba from the Telecommunications and Remote Sensing Laboratory, Pavia University, Italy³⁴). The number of spectral bands is 103, with image size of $610\times 340$ pixels. The spatial resolution is as high as 1.3 m per pixel, in contrast to the 20 m of the AVIRIS 92AV3C dataset. The dataset is also accompanied with a ground-truth map, which categorized the whole picture as nine types of materials, such as asphalt, meadows, gravel, trees, painted metal sheets, bare soil, bitumen, self-blocking bricks, and shadows. Figure 8 shows four examples of spectral bands for Pavia University scene and its ground truth.

Fig. 8

Illustration of Pavia University scene. (a) Four spectral bands. (b) Ground truth.

Following the similar approach to process AVIRIS 92AV3C dataset, we also tested our method on the Pavia University scene dataset. Figure 9 shows the spectral bands in two spectral windows, where Fig. 9(a) contains the bands in the spectral window found by the proposed method and Fig. 9(b) contains the bands from a random spectral window. It can be seen clearly that the bands in Fig. 9(a) are more clear and contain more discriminatory information than those in Fig. 9(b). This result coincides with our previous observation on the AVIRIS 92AV3C dataset (see results in Fig. 5).

Fig. 9

Spectral bands in (a) the spectral window found by the proposed method and (b) a randomly sampled spectral window.

By using the newly constructed reference map, we also assessed the change of classification accuracy as spectral bands are progressively selected. Table 1 lists the classification accuracies using SVMs as classifiers and the band subsets are selected from 5 to 85 bands, respectively (see the first row of Table 1). The second row of Table 1 lists the classification results using the reference map estimated from the suboptimal spectral window found by the proposed method. The third row contains the results using a randomly selected spectral window, and the fourth row lists the results using the ground-truth survey. It can be seen from Table 1 that when only five bands are selected for classification, the proposed method obtained 63.02% accuracy, compared to 60.92% (using ground-truth map) and 59.15% (using a random spectral window) of the two benchmarked methods. When about half of the bands (45 bands) are selected, the proposed method achieved 86.18% classification accuracy, only 1.5% lower than that using all the 103 bands. The results showed that the majority of the discriminatory information has already been included during band selection and confirmed the effectiveness of the proposed method.

Table 1

Comparison of classification accuracy, Pavia University scene.

Regarding the above results in Figs. 6 and 7 and Table 1, an interesting question may arise: why the estimated reference map generated a better result than that using the reference map accompanied with the dataset. This may be explained by the completeness of ground-truth labeling and the accuracy of the underlying probability estimation. First, in 92AV3C dataset, the ground truth is designated into 16 classes and they are not all mutually exclusive.^34,35 But PDF requires that states of a random variable are mutually exclusive. In this case, using the data in the ground-truth map to estimate the reference map’s PDF may produce significant errors (mismatch of the PDF’s definition). On the other hand, using the proposed method to estimate PDFs may avoid this problem because different real values are used to label different classes and they are mutually exclusive. Second, the accompanied reference map actually only labeled about half of the ground truth, i.e., the 16 classes of vegetations. The other half of pixels, e.g., highway, rail track, etc., are all categorized as background, presumably, because they correspond to uninteresting regions or were too difficult to label. This may put a strong argument to adopt the estimated reference map in the MI calculation because the estimated reference map differentiates different classes by different pixels’ value and all pixels can be fully utilized. Third, in the accompanied reference map, each class is labeled by an integer. For example, in AVIRIS 92AV3C, 16 classes of vegetations are labeled by 16 integers from 1 to 16, and the background is labeled by 0. Consequently, the estimated probability density will distribute in the integer domain. On the other hand, the estimated reference map is obtained by averaging several real spectral images, and the classes are differentiated by the real numbers characterized by their particular reflectance values. Since the MI is used to measure the utility of each band that was valued at the spectral response, the probability distribution estimated from the real reflectance values will match the MI calculation more than that using the artificially labeled values.