*3.5. Spectral Unmixing*

Finally, in order to evaluate the fractional abundances, a valid and unique method is needed for comparing the performance of the selected bands using the different methods studied. In this study, the least squares method has been employed to solve the inverse problem. According to [33], which studied the different methods of the estimation of the fractional abundances, the fully constrained least squares (FCLS) is introduced as a proper method in this regard. In order to accurately estimate the fractional abundances of the endmembers, two constraints—namely, the abundance sum-to-one constraint (ASC) and the abundance non-negativity constraint (ANC)—have been applied to the linear mixture model. In order to use the FCLS method, the constituents of the imaging scene should be fully known. This issue is considered using the information obtained from the ground truth map or the endmembers extracted from the image in the supervised and unsupervised manners, respectively.

#### **4. Results and Discussion**

In this section, the performance of the proposed ISI-PS algorithm is evaluated using the both simulated and real hyperspectral datasets. In the first subsection, several simulated hyperspectral images were employed, which have been produced using some spectra of the USGS spectral library and different scenarios. The constituent spectra of imagery and their fractional abundances were exactly known for these datasets. Then, the effect of selecting the optimal bands using the proposed ISI-PS algorithm to reduce the SV and endmembers' correlation was evaluated on the AVIRIS hyperspectral images. Finally, the results of the proposed method were compared with the results of the SZU method [13], which only dealt with the SV, and the results of the MTD method [18], which only tried to select the independent bands in an unsupervised manner.

#### *4.1. Simulated and Real Datasets Used*

## 4.1.1. Simulated Dataset

In this research, several hyperspectral images with a dimension of 100 × 100 pixels and various spatial patterns were simulated using some spectra from the USGS spectral library. These spectra were selected from different combinations of minerals of hydro-thermal alteration zones in geological applications (Table 1).


**Table 1.** The constituent spectra of the simulated images.

Several conditions were considered to simulate these images so that the resulting imagery reflected the real conditions as much as possible. In order to reconstruct the spatial patterns that appear in nature, the functions of the HYDRA software package [34] have been used, in which the two functions Legendre and Gaussian were employed to generate the fractional abundances. The Gaussian function could be performed using four different modes to generate the spatial patterns, and an example of each mode is illustrated in Figure 4.

**Figure 4.** An example of the fractional abundances that were generated using: (**a**) Spherical Gaussian fields; (**b**) Exponential Gaussian fields; (**c**) Rational Gaussian fields; and (**d**) Mattern Gaussian fields functions.

In this study, the effects of the endmembers' variability and the illumination fluctuation due to the topography on the spectrum of objects have been modeled according to [35]. The SV was simply characterized by spectral shape invariance [36]. In other words, while the spectral shapes of the endmembers were fairly consistent, their amplitudes varied considerably over the scene. Accordingly, the spectral variability of the *i*th endmember in each pixel can be modeled as Equation (3).

$$\mathbf{m}\_{i} = \psi\_{i}\mathbf{m}\_{i}^{\odot} + \eta\_{i\prime} \tag{3}$$

where **<sup>m</sup>**◦*i* is an endmember that is selected from the spectral library and *ψi* ≥ 0 is a random factor that affects the spectral amplitude of this endmember equally in all bands. *ψi* is the variable for each endmember in each class. This factor was generated by defining the range of the spectral variations (i.e., (1 ± *sd*)) and using the normal distribution. In the experiments, a standard deviation (sd) of 0.05 was adopted for the spectral variations. η*i* is a random noise with a zero-mean, which was considered to model those variations that are not modeled by *ψ<sup>i</sup>*.

By substituting Equation (3) in the linear mixture model (Equation 1), an equation was achieved to model the spectral variability for each pixel (Equation 4).

$$\mathbf{y} = \mathbf{M}\,\boldsymbol{\upmu}\,\boldsymbol{\upalpha} + \sum\_{i=1}^{p} a\_i \boldsymbol{\uppi}\_i + \mathbf{n}\_i \tag{4}$$

where ψ ≡ *diag*(*ψ*1, *ψ*2, ... , *ψp*) is a *p* × *p* diagonal matrix.

In order to apply these factors, a threshold was considered for the standard deviation of variations (sd). Then, the factors were randomly generated with a normal distribution in the range of 1 ± *sd*, which were multiplied by the amplitude of the original spectra. These factors were different for each pixel and properly model the effect of SV. The range of (0.05–0.15) was employed in experiments as the standard deviation of spectral variabilities.

The effect of the illumination fluctuation due to the topography was similar for all bands, and could therefore be considered as additive noise [35]. In order to model this effect, the spectral features matrix was considered to be fixed, and the fraction of each pixel was multiplied by a factor ( *γ*). These factors were randomly selected with a normal distribution in each pixel in the range of (0.95–1.05). In order to simulate the effects of the instrumental noises, Gaussian noise with the zero-mean and different ratios was added to the simulated scenes. Therefore,

$$\mathbf{y} = \mathbf{M}\boldsymbol{\gamma}\,\boldsymbol{\Psi}\,\mathbf{a} + \boldsymbol{\gamma}\sum\_{i=1}^{p} \boldsymbol{a}\_{i}\boldsymbol{\eta}\_{i} + \mathbf{n} = \mathbf{M}\,\mathbf{a} + \mathbf{v},\tag{5}$$

Equation (5) is the model that was considered in the experiments in this study, and still was linear. The illumination fluctuation, the SV, and the instrumental noise of images could be modeled using this equation. The flowchart of generating the simulated images is illustrated in Figure 5.

**Figure 5.** Flowchart of generating the simulated images.

## 4.1.2. LTRAS Dataset

The Russell Ranch Sustainable Agriculture Facility is a unique 300-acre facility near the UC Davis campus dedicated to investigating irrigated and dry-land agriculture in a Mediterranean climate. The goal of this facility—which is known as Long-Term Research in Agricultural Sustainability (LTRAS)—was used to investigate the impact of external factors such as crop rotation, farming systems (conventional, organic, and mixed), and the inputs of water, nitrogen, carbon, and other elements on agricultural sustainability [37].

Currently, the Century Experiment contains ten systems, which are two-year rotations and include corn/tomato, wheat/tomato, wheat/fallow, and wheat/legume rotations. Additionally, a perennial native grass system and a six-year alfalfa–corn–tomato rotation were initiated in 2012. The arrangemen<sup>t</sup> of each farm in this facility is illustrated in Figure 6a [38], along with the type of its irrigation and fertilization.

**Figure 6.** Long-Term Research in Agricultural Sustainability (LTRAS) dataset: (**a**) Planting schedule [38]; (**b**) Ground-truth map of the LTRAS farms.

The hyperspectral data used in this study were captured by the AVIRIS sensor on 3 August 2013, from low altitude with a 3.2-m ground pixel size. The region of LTRAS comprises 199 lines by 217 samples, and was extracted from the original image that had been radiometrically and geometrically corrected. Its true and pseudo-color composites are illustrated in Figure 7. The dimension of each farm is approximately 64 × 62 m, for which regarding the date of imaging and the planting schedule in Figure 6a, its classes have been extracted as Figure 6b. Plots 4-5 and 8-9—which had been planted with corn—were excluded from the ground truth map due to harvesting; as were Plots 5-4, which had not been planted with tomato as per the schedule. Because of the imaging date, the plots of cover crop had no covers, and were therefore considered as bare earth.

**Figure 7.** (**a**) True color composite; (**b**) False color composite; and (**c**) Color composite of the first three components of the maximum noise fraction (MNF) transformation.

The impacts of different irrigation systems and fertilization, as well as the effect of crop rotation were considered in this facility. Therefore, this dataset was proper to evaluate the SV phenomena in vegetables. The SV of pixels for three classes of wheat, tomato, and corn are illustrated in Figure 8.

**Figure 8.** The spectral variability of pixels for three classes of: (**a**) Wheat; (**b**) Tomato; and (**c**) Corn.

## 4.1.3. Salinas Dataset

This dataset was collected by the AVIRIS sensor on 9 October 1998, over Salinas Valley, Southern California, and is available as the at-sensor radiance unit. The scene comprises 512 lines by 217 samples, with 160 spectral bands (after discarding noisy and water absorption bands) in the wavelength range of 0.4–2.5 microns. Its nominal spectral and radiometric resolutions were 10 nanometers and 16 bits, respectively. This image was captured from low altitude with a 3.7-m ground pixel size. Its false color composite is illustrated in Figure 9a.

**Figure 9.** Salinas dataset: (**a**) False color composite of the AVIRIS image; (**b**) Ground truth map.

This dataset was collected from an agricultural region, and its ground truth map has been gathered into 15 classes, as illustrated in Figure 9b. It includes vegetables, bare soils, and vineyard fields with sub-categories as follows. The sub-categories of broccoli and green weeds were distinguished, with one having smaller and fewer weeds, while two had taller and more weeds, with both categories mostly covering the soil. The romaine lettuce sub-classes have been defined based on the planting week and their growth rates, which have different covers on the soil.

The soil was categorized into three sub-classes: the fallow rough plow class had recently been turned with larger clumps and appeared to have more moisture, while the fallow class was plowed soil with smaller clumps, and the fallow smooth class had even smaller clumps. The stubble class comprised bare soil and straw, and could also be considered a sub-class of the soil group. In the vineyard group, the untrained vineyard and the untrained grapes sub-classes were actually similar to each other. In the untrained vineyard sub-class, vine had been grown on wooden and plastic posts, and their canopies had nearly covered the soil. The situation of the selected classes at the time of imaging is shown in Figure 10.

**Figure 10.** Photograph of the selected classes in the region of the imaging.

#### 4.1.4. Indiana Indian Pines dataset

The third dataset used in the experiments was collected by the AVIRIS sensor over the Indian Pines Test Site in Northwestern Indiana in 1992. This image has a size of 145 × 145 pixels, and was acquired over a mixed agricultural/forest area, early in the growing season. The spatial resolution is approximately 20 m, and the radiometric resolution is 10 bits. The image comprised 220 spectral channels in the wavelength range from 0.4–2.5 micrometers, nominal spectral resolution of 10 nanometers. Bands 1–2, 100–114, 147–167, and 216–220 were removed from the dataset due to the noise and the water absorption phenomena, leaving a total of 177 radiance channels to be used in the experiments. This scene contained two-thirds agriculture and one-third forest or other natural perennial vegetation. For illustrative purposes, Figure 11a shows a false color composition of the AVIRIS Indian Pines scene, while Figure 11b shows the ground truth map available for the scene, with 16 classes. In our experiments, we considered a real situation in which most of the similar classes were included in the evaluations. Hence, 12 classes with an adequate number of labeled samples were selected for the experiments.

**Figure 11.** Indiana Indian Pines dataset: (**a**) False color composite of the AVIRIS image; (**b**) Ground truth map.

#### *4.2. Experiments on the Simulated Dataset*

The proposed ISI-PS algorithm was firstly performed on the simulated dataset, which was generated using the elements of Table 1. This data contained five datasets with different numbers and types of endmembers. In order to generate the fractional abundance for each dataset, a different pattern was employed according to Table 1. The main objective of this experiment was to evaluate the performance of the proposed method to deal with the endmembers' SV and decreasing the correlation of endmembers by selecting the optimal bands to generate accurate fractional abundances. Besides, the quality of endmembers' sets—which had been extracted in an unsupervised manner according to Figure 1—were compared with a spectral library that existed for the ground truth maps.

Several scenarios have been designed in order to evaluate the performance of the proposed method. In addition to the variety of spectral features and spatial patterns, different signal-to-noise ratios (i.e., 20:1, 25:1, and 30:1) have been employed to generate the simulated images. Equation (5) was used directly to simulate the first four datasets. In the case of the last dataset, in addition to Equation (5), two spectra from different species of one material were used to complicate the SV condition (Figure 12).

**Figure 12.** Endmembers used to simulate the last dataset.

As was mentioned in the fifth step of the proposed method, a threshold value (T) was employed to identify the independent bands in the PS. In other words, if the angle of the candidate band from all previously selected bands was greater than the threshold value, this band had a distinct behavior dealing with other selected bands.

The value of this threshold (T) was affected by the number, spectral similarity, and variety of endmembers in the scene. However, if a precise ground truth map were available, the root mean square error (RMSE) resulted from the comparison of this map, and the estimated fractional abundance could be used to properly evaluate this threshold value. In this regard, by increasing the threshold value in the range of 0.5–5 degrees with an increment of 0.25 degrees and evaluating the precision of the resulted fractional abundances, the threshold that led to the minimum error was selected. The results of different scenarios are provided in Table 2, along with the number of selected bands in each method.

If there were no in situ information for establishing a spectral library of the endmembers' variability, this information could be directly extracted from the image according to the early stages of the proposed algorithm. In this process—which was developed by applying some revisions to the SSPP algorithm—the virtual dimensionality (VD) of data was firstly estimated using the signal subspace identification algorithms (e.g., Hysime [28]) to be used for spectral dimension reduction by employing PCA or MNF transformations. Thereafter, the data used were also reduced in the spatial domain using spectral purity indices such as the PPI, and consequently, the number of candidate pure pixels was decreased. Using the mean spectrum of spectral clusters as the indicator of those classes without eliminating impure pixels led to the mean spectrum being affected by these pixels. In this case, the spectral changes among the clusters' pixels were not only due to the SV phenomena. By clustering those pixels that were probably pure, besides the tending of the class' mean towards the purity, the separability of classes would be properly shown.

Some other factors that challenged the performance of endmember extraction algorithms were sensitivity to noise, unusual pixels resulting from inaccurate atmospheric correction, and the image's hot spots. Therefore, extracting the homogeneous area of the image as a possible location for pure pixels could help to reduce the impact of these annoying phenomena. However, this could destroy the information of anomaly classes.

In order to identify the homogenous regions of the image used, a Gaussian filter with the functional form of *g*(*<sup>i</sup>* , *j* ) = 1 2*πσ*<sup>2</sup> *e*<sup>−</sup> *i* <sup>2</sup>+*<sup>j</sup>* 2 <sup>2</sup>*σ*<sup>2</sup> was firstly applied to each band of the hyperspectral image (Equation (6)). In this filter, the parameter σ controls the amount of spatial smoothing.

$$\mathbf{X}\_{b}^{F}(i,j) = \sum\_{i'=-a}^{a} \sum\_{j'=-a}^{a} \mathbf{g}(i',j') \cdot \mathbf{x}\_{b}(i+i',j+j') \tag{6}$$

where **X***Fb* (*i*, *j*) is the value of the pixel (*i*,*j*) in the band b of the filtered image. The value of *a* = (*w* − 1)/2 was determined regarding the filter dimension *gw*<sup>×</sup>*w*, and showed the spatial index of the filter. In order to generate the spatial homogeneity index (SHI), the RMSE of the discrepancies of the original image **X** and the filtered imaged **X***<sup>F</sup>* was computed using Equation (7) [27].

$$RMSE \left[ \mathbf{X}(i, j), \mathbf{X}^{\overline{F}}(i, j) \right] = \left( \frac{1}{B - 1} \sum\_{b = 1}^{B} \left( \mathbf{X}\_b(i, j) - \mathbf{X}\_b^{\overline{F}}(i, j) \right)^2 \right)^{\frac{1}{2}} \tag{7}$$

In this way, a layer was obtained where the value of the pixel (*i*,*j*) indicates the homogeneity of that region of the image. The lower the pixel value, the more homogeneous that region will be.

Finally, pixels that were located in each cluster were sorted based on the spectral purity index, and the purest ones that were located in homogenous areas were selected as the indicators of those clusters. In the experiments, 20 percent of the pixels of each cluster and 30 percent of the homogenous pixels with the best scores were employed. The results obtained from the subsections of this process are illustrated in Figure 13 for the first dataset.

**Figure 13.** The results obtained from the subsections of the process of extraction of the spectral library of each class from the image: (**a**) Output of the pixel purity index (PPI) for greater than zero pixels; (**b**) The result of band filtering; (**c**) The clustering of PPI's output; (**d**) Map of the homogeneity scores of pixels; (**e**) 20 percent of pixels of each cluster with the maximum value of PPI; (**f**) 30 percent of pixels of each cluster with the best homogeneity; (**g**) Fusion of maps (**<sup>e</sup>**,**f**) to extract the final pure pixels; (**h**) Ground truth map of the image's endmembers.

The results obtained from the SU of the five simulated datasets are provided in Table 3. It is worth mentioning that the endmembers' sets were directly extracted from the image in these experiments. For comparison purposes, the data used in these experiments were quite similar to those that were used in the previous supervised experiments.


Rational 5 30:1 70 44.58 0.61 0.071 109 62.38 **0.58** 0.074 23 4.00 **36.43** 0.60 **0.068** 224 76.40 0.62 0.075 Mattern 7 25:1 61 36.41 **0.19 0.043** 109 45.14 0.25 0.044 40 3.75 **33.31** 0.23 **0.043** 224 46.38 0.21 0.044 Mattern 10 30:1 55 43.32 0.55 0.175 152 70.40 0.59 0.175 12 3.00 **34.58 0.53 0.174** 224 74.47 0.60 0.176 #S is the number the image's components; #F is the number of selected bands; SNR is the signal-to-noise ratio; Cond is the Condition number; Corr is the mean Correlation of

endmembers; T is the angle threshold and bold number is the best result.

**Table 2.** The accuracy assessment results of the selected bands by maximum tangent discrimination (MTD), stable zone unmixing (SZU), and the proposed instability index (ISI)-prototype space (PS) algorithms, in comparison with using all bands to estimate the fractional abundances through a supervised spectral unmixing on the

132

According to Tables 2 and 3, the proposed ISI-PS algorithm always provided proper accuracy in the estimation of fractional abundances of endmembers by selecting the optimal bands. Moreover, the results obtained from the supervised and unsupervised experiments had good agreemen<sup>t</sup> with each other. Therefore, the spectral library of endmembers' variability was correctly extracted. The accuracies of the fractional abundances obtained from the MTD and ISI-PS methods were compatible. However, by comparing the position of selected bands in these two methods, it is obvious that in the MTD algorithm, the SV of endmembers was neglected, and the most separable bands were selected only by considering the spectral feature of each class. When the SV of the endmembers' sets was not to the extent that the spectrum of classes highly conflicted in the overlapping regions, the results of the two methods were close to each other. However, if the SV disrupted the separability of classes, the proposed ISI-PS method led to more accurate results by selecting the bands with the minimum spectral conflicts.

#### *4.3. Experiments on Real Datasets*

In order to evaluate the effect of SV on the results of unmixing, the ground truth map of the data used was needed, as well as the spectral library of the intra-class variations of each endmember. Providing the sub-pixel fractional abundance of the image's components was practically impossible, and collecting a spectral library from the variation of each component was an expensive and time-consuming process. However, in the case of the LTRAS dataset, for which only one crop was planted in each farm, the fraction of the components could be supposed to be 100 percent in the related plot; and observed spectral variations among the pixels of each plot could be seen as their SV due to the factors mentioned in Section 4.1.2.

However, it is worth mentioning that the change of the ratio of plants and background soil in each farm will lead to the change of the received spectra. This type of change was not considered as spectral variability. In this study, regarding the rather limited area and homogeneity of the farms, it was assumed that the mixture of materials was accrued with a constant ratio.

The extraction of the statistics of each class was performed using five percent of its pixels, which were selected randomly. This subset was used as the training dataset, and the remaining part was employed as the test dataset. Regarding the intra-class variations, the spectral features of the training pixels were employed for establishing a spectral library for each class. The mean spectrum of these sets was then used as the spectral indicator of the related classes.

As previously mentioned, the correlation of bands was not considered in the prioritization process of the SZU method. The selected bands using this method are illustrated by green bars in Figure 14. As can be seen, several bands were selected in the vicinity of each other in the event that the behavior of the endmembers was close together in these bands. In other words, if this redundancy could be reduced by a meaningful selection of the optimal bands in each region, more accurate and computationally-effective results could be achieved.

In the proposed ISI-PS algorithm, the angle of bands was used in the PS to deal with the bands' correlation. The selection of a proper threshold (T) to eliminate the correlated bands was affected by the number, spectral similarity, and variety of endmembers in the scene, as well as the number of resulting bands. If a ground truth map were available, the proper angle to eliminate the redundant bands could be estimated using the estimation accuracy of the ground truth map. In this experiment, regarding the homogeneity of farms, a ground truth map was generated using the map shown in Figure 6b. In this regard, the fraction of each endmember in the related class was considered as one, and for other classes, the fractions were considered to be zero. The fraction of classes in each pixel was then estimated involving the selected bands in each method, and their box plots are illustrated in Figure 15. In this case, the closer the resulting fraction to one and the less the standard deviation of fractions, the better the performance of the algorithm for dealing with the intra-class variations. It is worth mentioning that the most SV had occurred in the corn, wheat, and tomato classes. As can be seen, in comparison with the MTD method, the proposed ISI-PS algorithm tried to reduce the SV in

these classes by decreasing the median and the standard deviation of the estimated fractions. However, in the other classes, the performances of the studied methods were close to each other.

**Figure 14.** Representation of the selected bands using the MTD, SZU, and the proposed ISI-PS methods.

**Figure 15.** The box plot of the resulting fractions for each class using the selected bands by the studied method, using: (1) all bands; (2) MTD bands; (3) SZU bands; and (4) ISI-PS bands.

In this experiment, the threshold value of the correlated bands was considered as 1.25 degrees. Figure 16a shows the RMSE of the fractional abundance estimation of the LTRAS classes using the bands obtained from the MTD, the SZU, and the proposed ISI-PS algorithms, for threshold values of 0.25–5 degrees by an increment of 0.2 degrees. Index 6 was equivalent to the threshold value of 1.25 degrees, and caused the selection of 46 bands from the original 170 bands. The selected bands using the MTD, the SZU, and the ISI-PS algorithms are illustrated in Figure 14 using red, green, and dashed blue lines, respectively. It is worth mentioning that the proposed ISI-PS algorithm provided the most accurate fractions.

The singularity of the coefficient matrix, which was generated using the selected bands, was evaluated using: (1) the condition number of the endmembers matrix and (2) the average correlation of the endmembers' correlation matrix. By increasing the number of selected bands using each method, these two measures are illustrated in Figure 17, as well as for the full dimension of the data (i.e., all bands).

**Figure 16.** The RMSE of estimating: (**a**) The fractions; (**b**) The LTRAS image by means of the selected bands using the MTD, the SZU, and the ISI-PS algorithms for different threshold values.

**Figure 17.** Plot of: (**a**) The condition number of the endmembers matrix; (**b**) The average correlation of the endmembers' correlation matrix for the different band selection methods by increasing the number of the selected bands.

Finally, to evaluate the role of the selected bands on the results, the original image was reconstructed using the same set of endmembers and fractional abundances obtained from each method. In other words, by considering the accuracy of endmembers and applying the LMM, each pixel (*y*) of the original image could be approximated using *y*ˆ = ∑*pi*=<sup>1</sup> *αi***m***i*, where *p* is the number of endmembers, *αi* is the estimated fractions for each endmember, and **m***i* is the *i*th endmember. Accordingly, the original and the reconstructed images using the LMM could be considered as *I* = (*yk*)*Nk*=<sup>1</sup> and *R* = (*y*<sup>ˆ</sup>*k*)*Nk*=1, respectively, where *N* is the number of pixels. The reconstruction accuracy could be estimated using Equation (8). The average reconstruction error of each method is illustrated in Figure 16b, and the obtained results from the SU process are provided in Table 4.

$$RMSE\ (l, R) = \left(\frac{1}{N - 1} \cdot \sum\_{k=1}^{N} [y\_k - \hat{y}\_k]^2\right)^{\frac{1}{2}}\,. \tag{8}$$


**Table 4.** The accuracies obtained from the selected bands using the MTD, the SZU, the proposed ISI-PS methods, and the full dimension real hyperspectral datasets to estimate the fractional abundances.
