Translator Disclaimer
11 January 2017 Modeling linear and intimate mixtures of materials in hyperspectral imagery with single-scattering albedo and kernel approaches
Author Affiliations +
Linear mixtures of materials in a scene often occur because the resolution of a sensor is relatively coarse, resulting in pixels containing patches of different materials within them. This phenomenon causes nonoverlapping areal mixing and can be modeled by a linear mixture model. More complex phenomena, such as the multiple scattering in mixtures of vegetation, soils, granular, and microscopic materials within pixels can result in intimate mixing with varying degrees of nonlinear behavior. In such cases, a linear model is not sufficient. This study considers two approaches for unmixing pixels in a scene that may contain linear or intimate (nonlinear) mixtures. The first method is based on earlier studies that indicate nonlinear mixtures in reflectance space are approximately linear in albedo space. The method converts reflectance to single-scattering albedo according to Hapke theory and uses a constrained linear model on the computed albedo values. The second method is motivated by the same idea, but uses a kernel that seeks to capture the linear behavior of albedo in nonlinear mixtures of materials. This study compares the two approaches and pays particular attention to these dependencies. Both laboratory and airborne collections of hyperspectral imagery are used to validate the methods.



In a broad sense, spectral mixing phenomenology in multispectral and hyperspectral imagery can be treated in two ways, depending on how the materials in a scene are presumed to be mixed. A linear mixing model is suitable in cases where materials are presumed to be nonoverlapping areas and can be mathematically expressed as a weighted linear combination, where the weights are associated with the abundances of each material. The endmembers are spectra (real or ideal) representing unique materials in a given image such as water, soil, and vegetation. Abundances are the percentage of each endmember within a given pixel. On the other hand, it is well known that intimately mixed materials frequently exhibit nonlinear spectral mixing. Granular materials, such as soils, are often intimate mixtures of numerous different inorganic (minerals) and organic (humic) substances. And since soils are often significant constituents of spectral remote sensing scenes, intimate mixing may safely be assumed to be a common phenomenon. Another case is the mixing of fluids, such as oil and water, as occurred during the Deepwater Horizon oil spill disaster in 2010. We will consider this second case in more detail later.

Linear mixing is modeled as a linear combination of the spectra of multiple endmembers. The problem is posed in a number of ways using theory from linear statistical models with variations that impose either physical or sparseness constraints.17 Sites in a scene are labeled as containing one, two, or perhaps many endmembers. Some approaches begin with what is sometimes referred to as a full model, where the full model contains all possible variables (endmembers), and subsequently eliminates variables that do not contribute to the statistical significance (e.g., using an F-statistic) of the model.8 Alternatively, other approaches are synonymous with stepwise regression, where the process begins with a pair of variables (endmembers) and introduces new variables if they contribute significantly to the model.6,8

Recent research on linear unmixing has resulted in a number of alternative approaches. A stochastic mixture model (SMM) has been proposed with the goal to incorporate inherent endmember variance.9 Related to the SMM is a continuous version of the method, known as the normal compositional model (NCM).10 In the NCM, the mean vectors and covariance matrices of the endmembers are assumed adequate to capture their spectral variability. However, a big challenge to the NCM is parameter estimation. Very recently, a Bayesian algorithm called “normal endmember spectral unmixing” has been proposed to enhance parameter estimation performed in the NCM.11 Research has also been performed on the use of sparse modeling for modeling linear mixtures in hyperspectral imagery,12,13 a data-driven stochastic approach,14 as well as some advanced methods of nonnegative matrix factorization (NMF), such as projection-based NMF,15 and an adaptive L1/2 sparsity-constrained NMF.16

Intimate mixtures exhibit nonlinear spectral mixing behavior and are modeled as nonlinear combinations of spectra from multiple endmembers. In cases of intimate mixing, a linear spectral unmixing inversion applied to such nonlinear mixtures will yield subpixel material abundance estimates that do not equal the true values of the mixture constituents. An example of this is provided by Keshava and Mustard.17

An intimate mixture model can be described by nonlinear functions, which are justified by Hapke scattering theory18 and photometric phase functions.1921 In such an approach, reflectance is converted to single-scattering albedo (SSA). An exact model can be difficult to obtain, but can show up to 30% improvement in measurements over the linear mixing model when intimate mixtures are present.22 Many years ago studies showed success in using a constrained energy minimization method and other linear methods applied to reflectance spectra transformed to SSA data.2325

More recently, several alternative methods of spectral unmixing of nonlinear spectral mixtures have been proposed. Reviews of nonlinear spectral mixture analysis are given in Refs. 17, 26, 27, and 28 and references cited therein.

Examples include nonlinear unmixing based on sparse NMF,29 artificial neural networks with polytope decomposition,30 harmonic mixture modeling,31 and kernel-based methods.3239

Kernel functions provide a way to generalize linear algorithms to nonlinear data.32,33 In the cases of detection and classification applications, they can induce high-dimensional feature spaces. In these spaces, previously nonseparable classes are made linearly separable. Thus, linear methods can be applied in this new feature space that provides nonlinear boundaries back in the original data space. Another example is the kernel principal component analysis method.34 The kernel, in this case, is not used to induce a high-dimensional space, but is used to better match the data structure through nonlinear mappings. It is in this mode that kernels can be applied to obtain abundances in nonlinear mixtures while essentially using a linear mixture model. What is more appealing is that the physics suggests that such a method is ideal if one can model the kernel correctly.

A limitation of earlier kernel algorithms for material detection and scene classification is that they produced abundance estimates that did not meet nonnegativity and sum-to-one constraints. This was overcome by the development of a kernel fully constrained least squares (KFCLS) method that computes kernel-based abundance estimates to meet the physical (nonnegativity and sum-to-one) abundance constraints.35 Further investigation of the KFCLS method has resulted in (1) the development of a generalized kernel for linear and intimate (nonlinear) mixtures37 and (2) an adaptive kernel-based technique for mapping linear and intimate nonlinear mixtures.38 The generalized kernel and adaptive techniques provide a way to adaptively estimate a mixture model suitable to the degree of nonlinearity that may be occurring at each pixel in a scene. This is important because a scene may contain both linear and intimate mixtures and it is not always known a priori which model is appropriate on a pixel-by-pixel basis. This problem was investigated in a study performed by Broadwater and Banerjee.38 Building upon their work, a study investigating the behavior of the generalized KFCLS and adaptive kernel-based techniques was performed using both user-defined and automatically generated endmembers determined by the support vector data description algorithm.39

An investigation using laboratory data was also performed comparing the performance of the fully constrained least squares (FCLS) method applied to spectra converted to SSA with the generalized KFCLS applied to reflectance spectra and another recently published kernel-based method, the “K-Hype” kernel applied to reflectance spectra.40 One conclusion of this study: similar accuracy in abundance estimates can be achieved using the SSA-based method as compared to the generalized KFCLS, but with a much faster computation time. However, not all kernel methods are the same. The study also determined the “K-Hype” kernel method was not worth pursuing further for this type of analysis.

Overall, the impact of nonlinear spectral mixing on algorithm results must be well understood if we are to achieve a major goal of hyperspectral imaging (HSI): the geospatial mapping and quantification of materials that comprise remotely sensed scenes. In this study, we further this understanding by investigating both phenomenology-based SSA and mathematical-based kernel methods; and we investigate some issues noted while conducting earlier studies.39 The experiments include both laboratory and airborne HSI.

There has been a relative paucity of real-spectral image data sets that contain well-characterized intimate mixtures that can be used for algorithm development and testing. An important aspect of this study is the generation of a set of HSI data that are truly nonlinear spectral mixtures for which very precise truth information exists. Our laboratory experiment is performed on highly controlled data containing predetermined, nonlinear (intimate) mixtures of two materials. In addition, airborne experiments are performed on a scene of HyMap data collected over Oahu, Hawaii, and a scene collected by the airborne visible/infrared spectrometer (AVIRIS) instrument over the oil spill region in the Gulf of Mexico during the Deepwater Horizon oil incident.41,42


Description of Algorithms


Fully Constrained Least Squares and Linear Kernel-Based Mixing Models

The fully constrained least squares (FCLS)7 spectral mixing model can be expressed as

Eq. (1)

x=Ea+ϵwith the constraintsa(n)0,  nandn=1Na(n)=1,
where x is an L×1 vector containing the spectral signature of the current image pixel, E is an L×N matrix containing the endmember signatures thought to be within the scene (the n’th column contains the n’th endmember spectrum), and a is an N×1 vector containing the abundances [the n’th entry represents the abundance value a(n)]. ε is an L×1 vector containing the error. L is the number of bands in an image spectra and N is the number of endmembers used in the model. The abundance estimates a^ for the model can be computed through constrained least squares regression.

Alternatively, the linear model can also be posed as a special case of the kernel-based mixing model derived previously by Broadwater et al.,35 where the abundances of mixture components are estimated through the process

Eq. (2)

using the linear kernel K(x,y)=xty with the same constraints as Eq. (1). The abundance estimates are computed using a quadratic programming method to enforce the nonnegativity constraint.

Either way, once a^ is obtained, the estimator vector for a mixed spectral signature is

Eq. (3)

and the estimator for the error is

Eq. (4)


The root mean-squared error (RMSE) between the observed and estimated spectral signature is measured as

Eq. (5)


FCLS has been shown to be successful for modeling linear mixing phenomenology.7 In the present study, we will be using FCLS in two ways. The method will be used as a benchmark to compare with the proposed nonlinear methods. The method will also be used in one of the nonlinear approaches, where we will apply the FCLS to spectra that have been converted to SSA, as discussed below in Sec. 2.2.


Using Hapke Theory: Fully Constrained Least Squares Applied to Single-Scattering Albedo Spectra

Previous studies indicate that intimate (nonlinear or macroscopic) spectral mixtures in reflectance space may be linearized when transformed into SSA space. In order to understand this, we first mention the relationship between reflectance and SSA. Reflectance is the ratio of reflected radiance to incident irradiance.43 SSA is the ratio of the total amount of power scattered to the total power removed from the wave. If we assume the materials of interest are Lambertian, then we can consider the case of hemispherical reflectance and the mathematical relationship between reflectance and SSA can then be written

Eq. (6)

where r is the reflectance, H is Chandrasekhar’s function for isotropic scattering, w is the SSA, μ0 is the cosine of the angle of incidence, and μ is the cosine of the view angle.19

For a number of intimately mixed materials, the average SSA can be written as

Eq. (7)

where for material k, the SSA of the material is wk, the mass fraction is Mk, the diameter of the particles of material is dk, and the solid density of the particles is ρk.38

According to Eq. (7), intimate mixtures are a linear mixture of the albedos wk and these are nonlinearly mapped to reflectance by Eq. (6). This linear mixing of the SSAs can be seen trivially, as we can write Eq. (7) as

Eq. (8)


Eq. (9)

Fk is the relative geometric cross section. 38 Noting the form of Eqs. (8) and (9) as compared to Eq. (1) for linear mixing, the terminology for relative geometric cross section and “abundance” can be used interchangeably.

Accordingly, we investigate this behavior by applying a linear mixing method on albedo; specifically, by applying the FCLS method on data that has been converted to SSA.

Conversion to SSA is described further in Resmini et al.23 and Resmini.24 Both studies follow Hapke;18 and Mustard and Pieters1921 assuming the reflectance spectra are bidirectional. The expressions to transform reflectance spectra to SSA are given by Eqs. (10) and (11) for bidirectional (bd) reflectance and for hemispherical-directional (hd) reflectance, respectively. In the derivation of both expressions, the phase angle is large enough that the opposition effect is assumed negligible.

Eq. (10)


Eq. (11)

In Eqs. (10) and (11), Γ is the reflectance factor (see Sec. 2.1, Hapke18); r, w, μ0, and μ are the same as in Eq. (6). The two different equations can be used to generate the two sets of SSA spectra for bd and hd reflectance.

In the experiments that follow, we refer to this approach as the “FCLS on SSA” method, or simply denote it as the SSA method. The conversion of reflectance spectra to SSA using Eqs. (10) and (11) is very fast as compared to the generalized kernel least squares (GKLS) method.


Generalized Kernel Fully Constrained Least Squares (Fixed and Automated)

Choosing the linear function K(x,y)=xty for the kernel in Eq. (2) works well for modeling linear mixtures; however, it is not a suitable kernel for intimate mixtures. A physics-inspired kernel was proposed and shown to provide significantly improved behavior to model nonlinear mixtures, but a result of that effort was that although each kernel provides good results for the type of mixing intended, only one kernel or the other could be used, for either linear mixtures or intimate mixtures, but not both.36

The kernel approach was further developed by Broadwater and Banerjee37,38 into a generalized method for adaptive linear and intimate mixtures. This method is motivated by attempting to simulate Hapke theory for SSA by making use of the kernel

Eq. (12)

The kernel in Eq. (12) can model either linear or intimate mixtures as long as the correct value of γ is used. If γ is very small then Kγ(x,y) approximates linear mixing. If γ is large then Kγ(x,y) approximates intimate mixing in cases when the reflectance occurring from intimate mixing is modeled as in Eq. (6).

In this study, we investigate a fixed-γ generalized kernel least square” (fixed-γ GKLS). The computation is similar in form to Eq. (2) except the minimization is done according to

Eq. (13)

using Eq. (12) for a user-defined γ, and the same constraints as Eq. (1).

For a specified γ, we implement this by transforming the observed x and E into kernel space

Eq. (14)


Eq. (15)

where xγ and Eγ are presumed to satisfy the linear model described by

Eq. (16)

xγ=Eγa+  ϵγ.

Equation (16) has the same constraints as Eq. (1) and the abundance estimates a^ can be computed in the same manner as Eq. (1). The estimated mixed spectra in kernel space is computed as

Eq. (17)

x^γ=  Eγa^.

The estimated error in kernel space is

Eq. (18)


The estimated error vector ϵ^ in the original space is computed as Eq. (4), where the components of the estimated mixed spectra x^(n) are computed using the inverse transform of Eq. (14). The RMSE is computed by Eq. (5).

In addition to a fixed-γ GKLS implementation, an automated GKLS method is investigated that attempts to select the most appropriate gamma according to

Eq. (19)

where a^γ is the abundance estimate and Kγ(x,y) is the kernel evaluated with the parameter value γ. A numerical optimization based on the golden search method can be used to minimize Eq. (19).44 In computing the optimization, we seek to achieve a minimum of the model’s RMSE. In this manner, at least theoretically, the automated GKLS method has the ability to respond differently to differing degrees of nonlinearity to select the best gamma and compute more precise estimates of abundance.


Description of Experiments


Overview of Experiments

Three experiments are performed in this study. The first is an experiment conducted in the laboratory. This is an important aspect of our investigation: a set of HSI data is generated that are truly nonlinear spectral mixtures for which very precise truth information is documented. The laboratory experiment is performed on highly controlled data containing predetermined, nonlinear (intimate) mixtures of two materials. In this experiment, two granular materials were custom fabricated and mechanically mixed to form intimate mixtures. The materials are spherical beads of didymium glass and soda-lime glass, both ranging in particle size from 63 to 125  μm. The glass beads materials are both translucent. Their chemical composition, densities, and particle size range are well known. The mixtures, which exhibit largely nonlinear spectral mixing, were then observed with a visible/near-infrared (VNIR; 400 to 900 nm) HSI microscope. The second and third experiments are conducted in the field with airborne hyperspectral imagery: the second experiment is performed using data from the HyMap sensor acquired during a campaign by the Navy Research Laboratory (NRL) over Oahu, Hawaii, from January 24 to February 1, 2009. The third experiment utilizes AVIRIS data collected over the oil spill region in the Gulf of Mexico during the Deepwater Horizon oil incident. The SSA, FCLS, and GKLS methods are each applied to data from the three collections.


Experiment 1: Laboratory Experiment

A Resonon Pika II imaging spectrometer with a Xenoplan 1.4/23-0902 objective lens shown in Fig. 1(a) is used to measure the didymium and soda-lime glass bead mixtures.45 We also used an Edmund Optics Gold Series 1.0× telecentric lens. However, this lens gives 8  μm/pixel data, providing a spatial resolution higher than what was required for the analyses. Consequently, the Edmund data are not reported in our results.

Fig. 1

The laboratory apparatus is shown: (a) a photograph of the VNIR HSI microscope. The Resonon Pika II is used with the Xenoplan 1.4/23-0902 objective lens; (b) 35-mm digital single-lens reflex camera photograph of the 96-well plate containing the glass beads. (Source of photo: taken by coauthor Dr. David Allen, NIST, property of U.S. government).


The Pika II is mounted nadir-looking at a mechanical translation table on which the sample to be imaged is placed. This device is a pushbroom sensor with a slit aperture, thus the need for a translation table to move the sample to facilitate hyperspectral image cube formation. The height of the sensor above the table is user selectable; a height was chosen such that all mixtures are captured in the same scene so that the data have a ground sample distance of 75  μm/pixel. Though capable of acquiring 240 bands from 400 to 900 nm, the sensor was configured to acquire 80 bands by binning (spectrally by three) resulting in a sampling interval of 6.25  nm and high signal-to-noise ratio spectra. Four quartz–tungsten–halogen (QTH) lamps are used for illumination approximating a hemispherical-directional illumination/viewing geometry.

Sensor and translation table operation, data acquisition, and data calibration are achieved by software that runs on a laptop computer. Calibration consists of a measurement of dark frame data (i.e., acquiring a cube with the lens cap on) and a measurement of a polytetrafluoroethylene (PTFE) reference plaque (large enough to entirely fill the field-of-view). Then, for each HSI cube measured, the sensor’s software first subtracts the dark data and then uses the PTFE data (also dark subtracted) to ratio the spectral measurements to give relative reflectance (also known as reflectance factor: Hapke18 and Schott43).

Three binary mixtures (and the two endmembers) are constructed and emplaced in the wells of a 96-well sample plate: 0/100%, 25/75%, 50/50%, 80/20%, and 100/0% of didymium/soda-lime (percentages by volume). This was done as follows: five cells of a 96-well sample plate, spray-painted flat black, were filled with the various glass bead mixtures; this is shown in Fig. 2. The volume of each cell is 330  μL (0.33 mL). Three binary mixtures and the two 100% endmembers are constructed and emplaced in five of the wells of the 96-well sample plate: 0/100%, 25/75%, 50/50%, 80/20%, and 100/0% didymium/soda-lime. The glass beads and their mixtures also display subtle, though interesting, gonioapparent changes in color. We remark that the glass bead particle size range is much larger than the VNIR wavelengths used in this analysis.

Fig. 2

RGB composite images of the hypercube showing the training and test regions used in the experiment trials. Figure 2(a) is an RGB composite of the HSI cube showing training regions. Figure 2(b) is an RGB normal color composite of the HSI cube showing test regions.


Fifteen data Pika II HSI cubes were acquired; however, we focus here on the analysis of one cube comprised of 640 samples, 400 lines, and 75 bands ranging from 0.434 to 0.885  μm. Of the 80 bands acquired the first five were discarded due to noise content.

Training and test data were extracted from the selected hyperspectral cube. Figure 2 shows polygons defining the training and test regions drawn on top of a red–green–blue (RGB) color composite context image. The training regions are shown in Fig. 2(a) and the test regions are shown in Fig. 2(b). Note that none of the test regions overlapped the training regions. Figure 2(a) shows the small training polygon regions defining the three training endmembers: DiDy (100%), lime (100%), and background. Figure 2(b) shows the test areas used to quantitatively measure the performance of the algorithms. Five regions were defined, corresponding to the five mixtures 100% DiDy, 75/25% DiDy/lime, 50/50% DiDy/lime, 25/75% DiDy/lime, and 0/100% DiDy/lime. Two additional test regions of background spectra were also extracted.

Figure 3 shows plots of the mean vectors of the spectral data extracted within the training polygons shown in Fig. 2(a). These image-derived training endmembers, as just described, are used to investigate the three methods described in Sec. 2: (1) FCLS applied in reflectance space, (2) GKLS applied in reflectance space, and (3) FCLS applied in SSA space. For purposes of conciseness in reporting results, these methods henceforth will be referred to as the FCLS, GKLS, and SSA methods, respectively.

Fig. 3

Mean spectra for the three endmembers (DiDy, lime, and dark background) are shown.


Numerous factors affect the performance of the methods. In particular, three factors affecting the performance of all the methods are: (1) the number the endmembers used in a model, (2) how well these endmembers span the space in which the mixing occurs, and (3) the RMSE threshold for eliminating bad fits between the observed and model-estimated spectra. In addition, the GKLS method uses a parameter γ to determine the nonlinear behavior introduced by the kernel. We test the GKLS method at fixed values of γ: 0.1, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, and 6.0, as well as the automated GKLS method. For the SSA method, performance may be affected by the type of SSA conversion: hemispherical or bidirectional. If bidirectional, the view and illumination angles are factors. Although the laboratory instrument has approximately hemispherical illumination, we will conduct trials for SSA conversions made assuming both bidirectional reflectance with nadir view and illumination angles as defined by Eq. (10), denoted as SSA, as well as hemispherical reflectance as defined by Eq. (11), which we will refer to as SSA-H.

The algorithms are applied to the entire scenes. Both qualitative results (shown by images of the algorithm output) and quantitative results (applied in the test regions) are given.


Experiment 2: Airborne Experiment in Hawaii

Experiment 2 tests the proposed methods using a cube of HyMap data collected during a campaign by the NRL over Oahu, Hawaii. Flight line Az148 (acquired on February 1, 2009) was selected and analysis was performed on a subscene of the original scene with a size of 400 samples by 250 lines over government-owned property near the Waimanalo region of Oahu. The imagery consists of 125 bands, has a spectral range of 0.45 to 2.48  μm with a spectral resolution of about 0.01  μm, and a ground sampling distance (GSD) of 4  m. The calibrated, at-aperture radiance data were converted to reflectance by the NRL using their Tafkaa,46 atmospheric compensation tool. Tafkaa is a heavily modified version of ATREM.47 The reflectance spectra are scaled by a factor of 10,000 and converted to 16-bit integers. An RGB normal-color composite of this scene is shown in Fig. 4.

Fig. 4

An RGB composite of the HyMap scene used in the experiment. This subset contains 400 samples by 250 lines. [Source of imagery: U.S. Naval Research Laboratory (NRL).]


Ten fabric panels of varying size were placed at three different sites (and in varying local backgrounds) prior to the image collection. These panels are impermeable Tyvek-like house wrap, near white in appearance. The sizes of these panels are listed in Table 1. Figure 5 shows aerial photographs of the three sites. The locations of the panels at each site are annotated. Numerous ground photos were also taken of the panels. Figure 6 shows ground-truth photos for three of these panels (TY-4, TY-5, and TY-7).

Table 1

Target sizes: the sizes of the 10 target panels are listed in this table. The measurements were made in units of feet, but converted to meters for easier comparison to the estimated 4.0 GSD of the airborne HyMap scene.

TY012.7  m×3.0  mTY065.5  m×5.5  m
TY022.7  m×1.4  mTY072.7  m×3.2  m
TY032.7  m×1.5  mTY082.7  m×3.4  m
TY042.7  m×3.2  mTY092.7  m×7.6  m
TY052.7  m×3.0  mTY102.7  m×7.6  m

Fig. 5

Aerial photos of the panel study sites are shown: (a) the configuration of six Tyvek panels at site 1, denoted TY-4, TY-5, TY-6, TY-7, TY8, and TY-9; (b) the three Tyvek panels at site 2 denoted TY-3, TY-2, TY-1; and (c) the Tyvek panel at site 3 denoted as TY-10. [Source of imagery: U.S. Naval Research Lab.]


Fig. 6

Ground photos of three Tyvek panels located at site 1: (a) TY-4, (b) TY-5, and (c) TY-7. (Source of photos: USGS and NGA.)


Site 1 contains 6 of the 10 panels (TY-4, TY-5, TY-6, TY-7, TY-8, and TY-9). TY-7 is expected to be difficult to detect because of its subpixel size in the HyMap scene and because it is partially occluded by trees. Notice that TY-7 is somewhat difficult to see in the aerial photo shown in Fig. 1; however, it is considerably more difficult to see in the HyMap data because of the reduced spatial resolution. TY-5 is expected to be difficult because of size and shading. Site 2 contains three Tyvek panels (TY-1, TY-2, and TY-3). Site 3 contains one Tyvek panel. During the campaign, spectral measurements were acquired by the USGS and NGA using an ASD spectrometer. However, these field measurements were not used for training in this study. Instead image-derived spectra were used.

Figure 7 shows four of the eight endmember spectra used in experiment 2. Similar to the procedure used in the previous experiment, training regions for eight classes of materials in the scene were defined for use as endmembers. Seven of these belonged to naturally occurring classes, such as sand, gravel, water, and different types of vegetation having regional extent (e.g., multiple pixels). These endmembers were computed as mean spectra within regions defined by polygon overlays (similar to the overlays shown in Fig. 2). However, as shown in Table 1, the Tyvek panels are comparably small and mostly subpixel in extent. TY-6 is the largest panel and is just slightly larger than a full pixel (5.5×5.5  m). The spectral vector from the center point of TY-6 was extracted and used as the Tyvek endmember.

Fig. 7

Plots of four (out of the eight) endmember spectra used in experiment 2 are shown. The Tyvek endmember is the spectrum extracted from the image at the center point of TY-6. The others are mean spectra computed over polygon regions extracted from the image.



Experiment 3: Airborne Experiment: Deep Water Horizon

Experiment 3 uses a scene that was extracted from AVIRIS data collected during a JPL/USGS airborne campaign during the Deep Water Horizon (DWH) oil spill incident.41,42 As a part of these studies, during the DWH campaign, USGS collected numerous oil samples for laboratory measurements.41 Although precise ground truth is not available due to the dynamically changing environment of an oil spill on water, previous studies using feature-based spectroscopy methods determined the locations of three predominant oil thicknesses in the region.42

Out of the numerous flight lines of AVIRIS acquired on May 17, 2010, run 11 of AVIRIS was selected, which includes the actual incident site of the oil spill. This run was collected at an altitude of 28,000 feet resulting in a GSD of 9  m. The AVIRIS scene is 677 samples by 16,835 lines with 224 bands of reflectance data from 0.365 to 2.49  μm. As part of the calibration to reflectance, the USGS scaled the imagery by a factor of 20,000. The experiment focused on a subset scene of size 500 pixels by 500 lines by 224 bands. Out of the 224 bands, 97 bands were excluded: The first 39 bands (0.365 to 0.714  μm) were excluded to focus the algorithm on the longer wavelengths away from the visible, the last 10 bands (2.407 to 2.496  μm) were excluded because of low signal-to-noise, and the remaining bands were excluded because of atmospheric absorption.

Figure 8 shows a true-color RGB composite of the subset scene. The scene is mostly water with large portions containing a great diversity of oil thicknesses. Different states of the oil can be seen visibly, appearing in different colors ranging from almost black, to dark brown, to light brown, and to orange.

Fig. 8

An RGB composite scene (500 samples by 500 lines) from a portion of AVIRIS Run 11 is shown in (a). The corresponding center wavelengths for the red (R), green (G), and blue (B) layers are 0.658, 0.541, and 0.472  μm, respectively. (Source of imagery: U.S. Geological Survey.)


Figure 9 shows laboratory-measured reflectance spectra of four oil samples at differing thicknesses. Absorption features due to C-H with band centers (bc) at 1.2 and 1.73  μm are labeled as F1 and F2.

Fig. 9

Spectra of three oil samples at differing thicknesses are shown. Absorption features at F1 and F2 are due to C-H with band centers (bc) at 1.2 and 1.73  μm.


For purposes of training and defining the endmembers, laboratory spectra for pure oil samples with 4-mm thickness and 1.85-mm thickness are used, as well as an oil/water emulsion having a thickness of 1.0 mm. Because much of the emulsion behavior is difficult to model in the laboratory, user-defined image spectra are also used: orange-colored spectra (in the RGB true-color composite) and water spectra with trace oil.

The FCLS, GKLS, and SSA methods are applied to the AVIRIS data, reporting on the results using the fixed-γ implementation and the SSA conversion made assuming bidirectional reflectance with nadir input and output angles.




Overview of Results

The results for the three experiments are presented below, showing tables with model diagnostics that provide estimated abundance values for the estimated models and RMSE values for the model errors; abundance maps that provide a visual representation of the abundances of materials found in the scenes; and graphs showing a comparison of the estimated (simulated) spectra estimated by the computed models and the observed spectra at selected points.


Results of Experiment 1

The results for experiment 1 (the laboratory experiment) are shown in Figs. 10 and 11, as well as Tables 1 and 2. Figure 10 shows RGB color abundance maps for four of the trials (red=DiDy, green=lime, blue=background). Qualitatively, Figs. 10(a) and 10(b) show poor correspondence to the known mixtures shown in Fig. 2. Figures 10(c) and 10(d) show much better correspondence to the actual abundances given in Fig. 2. The variations in color within the discs containing the three mixtures (75/25, 50/50, and 25/75) are noteworthy. These variations indicate the methods are detecting notable variance in the abundance proportions. This indicates that the attempt to prepare mixture proportions that are as uniform as possible were not completely successful. These variations are real and the abundance estimation methods (particularly, GKLS and SSA) are responding correctly. Static clinging and other interparticle interactions likely account for the clumping and variations observed. Nonetheless, overall population averages captured by regions of interest delineated in Fig. 2(b) covering a large portion of the individual cells yield averages close to those intended in the experiment design phase.

Fig. 10

RGB composite images showing color abundance maps for four of the trials (red=DiDy, green=lime, blue=background). Abundance maps are shown for (a) FCLS method (linear), (b) GKLS at γ=0.1 (linear), (c) GKLS at γ=5.0 (nonlinear), and (d) SSA method (nonlinear). Figures (a) and (b) show poor correspondence to the known mixtures shown in Fig. 2. Figures (c) and (d) show much better correspondence to the actual abundance values.


Fig. 11

The observed and estimated mixture spectra (averaged) of the 50/50% region using the FCLS, GKLS, and SSA-H methods. The y-axis of (a) and (b) is in reflectance units with a range 0.0 to 1.0; the y-axis of (c) is in albedo units with a range 0.0 to 1.0. The GKLS results are shown for γ=5.


Table 2

Abundance results: the average estimated abundances are listed for the FCLS, GKLS, and SSA methods in the five test regions. The “truth” for these regions (actual physically measured proportions by volume) is 1.0, 0.788, 0.505, 0.242, and 0.0 for DiDy100, DiD075, DiDy050, DiDy025, and DiDy000, respectively. Results for the GKLS method using a fixed gamma of γ=0.1, γ=0.5, γ=1.0, γ=2.0, γ=3.0, γ=4.0, γ=5.0, and γ=6.0 are given, as well as for the automated GKLS.


Table 2 lists the average estimated abundances for the FCLS, GKLS, SSA, and SSA-H methods in the five test regions. The “truth” (actual physically measured) percent by volume of DiDy for these regions varied slightly from the goal of 100%, 75%, 50% 25%, and 0%. In reality, these were measured as 100%, 78.8%, 50.5%, 24.2%, and 0% for DiDy100, DiD075, DiDy050, DiDy025, and DiDy000, respectively. Results for the GKLS method using a fixed gamma of γ=0.1, γ=0.5, γ=1.0, γ=2.0, γ=3.0, γ=4.0, γ=5.0, and γ=6.0 are given, as well as for the automated GKLS. This table shows FCLS to be poor at predicting the abundances for DiDy075 (93.11% predicted versus 78.8% truth) and DiDy050 (74.97% predicted versus 50.5% truth), as well as very poor at predicting DiDy025 (63.88% predicted versus 24.2% truth). The results of GKLS for small values of γ agree with theoretical expectations of an approximately linear model. Specifically, the prediction of GKLS at γ=0.1 is almost exactly the same as the FCLS method.

Out of the eight γ values tested, GKLS at γ=5.0 provides the closest prediction for DiDy50 (49.08% versus 50.5%) and is only slightly worse at predicting the correct abundance than GKLS at γ=6.0.

Figure 11 shows the observed and estimated mixture (averaged) spectra of the 50/50% region using the FCLS, GKLS, and SSA-H methods. Visually, we can see both the GKLS (γ=5) and SSA-H methods provide a better fit as compared to the FCLS method.

Table 3 lists the RMSE of the fit between the estimated and observed spectral mixtures for selected points in the scene. Except for DiDy at 25% (0.242), the RMSE errors for FCLS were not considerably larger than the errors for the other methods. Yet we know FCLS is poorly predicting the known abundances for these samples. We conclude RMSE is not necessarily a good indicator of a method’s accuracy to predict abundance.

Table 3

Model diagnostics: the RMSE results of the fit between the estimated and observed spectral mixtures are listed for selected points in the scene. The first column lists the location and planned percentage mix of DiDy for these points. The actual physically measured mixes were 100%, 78.8%, 50.5%, 24.2%, and 0.0%. Bold values indicate lowest RMSE value retrieved for the indicated pixel.

(585, 248) 100%0.01410.02930.03690.03150.0320
(400, 228) 75%0.02450.02440.02480.03510.0388
(418,239) 75%0.01830.01770.03240.02620.0283
(405,251) 75%0.01790.01960.02050.02550.0273
(324, 049) 50%0.03570.02250.02220.01290.0190
(313, 062) 50%0.03590.03180.03900.01180.0175
(328, 067) 50%0.03330.02940.05330.02230.0294
(223, 314) 25%0.05260.03800.03790.01400.0241
(224, 327) 25%0.04950.03830.03840.01860.0298
(246, 330) 25%0.04200.02890.06360.01210.0197

The results in Tables 2 and 3 also show that the automated implementation of the GKLS method was not as successful as the fixed-γ GKLS (γ=3, 4, 5, or 6) for estimating the correct abundance. The automated GKLS method attempts to select the most appropriate γ based on achieving a minimum of the model’s RMSE. This metric seems to respond to linear and nonlinear mixtures; however, it does not appear to be a reliable metric to determine the degree of nonlinear behavior. RMSE also cannot be used to achieve the most accurate estimate of abundance. Consequently, an alternative approach should be sought for implementing an automated GKLS.

Determining which value of γ gives the best result for this experiment is difficult. In terms of abundance, γ=5.0 clearly provides the best estimate for the DiDy50 mix (49% abundance). However, γ=6.0 provides slightly better estimates of abundance for the DiDy75 and DiDy25 mixes, (84% and 0.31%, respectively). In terms of RMSE, γ=5.0 provides a better fit than γ=6.0 in most (but not all) cases. Noting that γ=5 provides the highly close prediction of abundance for the DiDy50 mixture and the RMSE fits are better, we henceforth consider γ=5.0 to provide the best GKLS result as compared to the others (although only slightly better than γ=6.0).

The processing times on the data consisting of 256,000 pixel vectors (640 samples by 400 lines) were recorded for each of the methods. The FCLS, SSA, and fixed-γ GKLS runs were completed in 9 s, 9 s, and 12 s, respectively. The automated GKLS run finished in 228 s. Note the automated GKLS method was slower by a factor of 19.


Results of Experiment 2

Results for experiment 2 are shown in Tables 4Table 5Table 67 and Figs. 12 and 13, showing the behavior for the FCLS, GKLS, and SSA methods, quantitatively and visually. The SSA method in this experiment is synonymous to the SSA method discussed in the first experiment, which is defined by Eq. (10) assuming bidirectional reflectance with nadir view and illumination angles.

Table 4

Model diagnostics of panels for experiment 2. The RMSE values of the FCLS, SSA, and GKLS methods are listed for selected panel pixels in the scene, as well as the corresponding γ value chosen by the GKLS method at these points.

T01(369, 126)0.0070.0170.0076.440T06A(211, 62)0.0160.0180.0116.440
T02A(365, 128)0.0200.0220.0174.569T06B(211, 63)0.0000.0000.0005.126
T02B(365, 129)0.0210.0290.0203.413T06C(212, 63)0.0100.0170.0096.440
T03A(363, 182)0.0200.0140.0206.440T07(187, 99)0.0050.0180.0050.002
T03B(362, 182)0.0390.0240.0396.440T08(201, 64)0.0040.0140.0040.006
T03C(363, 181)0.0040.0040.0040.005T09A(163, 33)0.0070.0100.0072.838
T04A(209, 79)0.0050.0210.0050.002T09B(164, 33)0.0070.0240.0070.002
T04B(209, 80)0.0030.0070.0031.778T10A(43, 106)0.0080.0140.0083.459
T05A(212, 76)0.0080.0170.0075.321T10B(43, 107)0.0070.0100.0062.479
T05B(212, 77)0.0030.0140.0031.714T10C(44, 107)0.0060.0130.0060.002

Table 5

Model diagnostics of vegetation mixes for experiment 2. The RMSE values of the FCLS, SSA, and GKLS methods are listed for selected pixels containing vegetation mixtures in the scene, as well as the corresponding γ value chosen by the GKLS method at these points.

(277, 50)0.0030.0130.0032.272
(277, 60)0.0040.0120.0040.149
(248, 68)0.0060.0130.0062.415
(99, 79)0.0090.0180.0076.440
(134, 117)0.0110.0270.0114.810
(233, 141)0.0050.0140.0041.611
(234, 141)0.0060.0130.0052.685
(270, 178)0.0050.0130.0052.129
(274, 179)0.0060.0190.0061.011
(275, 180)0.0050.0120.0053.069

Table 6

Abundance results for Tyvek panels in experiment 2. The abundance values computed by the FCLS, SSA, and GKLS methods are listed for selected locations of the Tyvek panels.

Panel nameLocation (x,y)WaterDrkGravelSandPitGrass healthyTreesVegNrTyvekTyvekPanel nameLocation (x, y)WaterDrkGravelSandPitGrass healthyTreesVegNrTyvekTyvek
FCLST01(369, 126)0.040.650., 99)
SSAT01(369, 126), 99)
GKLST01(369, 126), 99)
FCLST03C(363, 181), 64)
SSAT03C(363, 181), 64)
GKLST03C(363, 181), 64)
FCLST04B(209, 80)0.000.520., 33)
SSAT04B(209, 80)0.010.360., 33)
GKLST04B(209, 80)0.000.470., 33)
FCLST05B(212, 77)0.530., 106)0.000.670.
SSAT05B(212, 77)0.320., 106)
GKLST05B(212, 77)0.490., 106)0.000.570.
FCLST06C(212, 63)0.330., 107)0.000.570.
SSAT06C(212, 63), 107)
GKLST06C(212, 63), 107)0.000.440.
The bold values emphasize the predicted abundance for “Tyvek,” the primary material of interest.

Table 7

Abundance results for vegetation mixtures in experiment 2. The abundance values computed by the FCLS, SSA, and GKLS methods are listed for selected locations of vegetation mixes.

Vegetation mix location (x,y)WaterDrkGravelGrass healthyTreesVegNrTyvekBlg01Vegetation mix location (x,y)WaterDrkGravelGrass healthyTreesVegNrTyvekBlg01
FCLS(277, 50), 141)
SSA(277, 50), 141)
GKLS(277, 50), 141)
FCLS(277, 60), 141)
SSA(277, 60), 141)
GKLS(277, 60), 141)
FCLS(248, 68), 178)
SSA(248, 68), 178)
GKLS(248, 68), 178)
FCLS(99, 79), 179)
SSA(99, 79), 179)
GKLS(99, 79), 179)
The bold values emphasize the predicted abundance for three types of vegetation, the primary materials of interest.

Fig. 12

A true color composite of the subset of the scene containing sites 1–3 is shown in (a) (source of imagery: NRL). The abundance maps corresponding to Tyvek are displayed in (b) and (c), respectively, for the GKLS and SSA methods. If the RMSE was greater than the chosen threshold value, RMSE_MAX=0.015, then the computed abundances for that pixel were rejected and set to zero. The circles image shows the location of the tyvek panels. A number of false alarms are also made at this RMSE threshold.


Fig. 13

A false-color RGB composite abundance map of the subset scene for the vegetation endmembers using the FCLS method is shown in (a). Composite RGB abundance maps using the GKLS and SSA methods are shown in (b) and (c), respectively. Notable differences in pattern between the GKLS and SSA methods can be observed. These maps are generated with the abundance layers corresponding to the grass (red), trees (green), and brush (blue). If the RMSE was greater than the threshold value, RMSE_MAX=0.015, then the computed abundances for that pixel were rejected and set to zero.


The abundance maps for experiment 2 are shown in Figs. 12 and 13. Figure 12 shows maps for the Tyvek panels using the GKLS and SSA methods. The map for the FCLS method is not shown because this map is almost identical to the GKLS method. The colored circles show where the panels are located. Due to limited dynamic range in the printed version, some of the low abundance panels (T02, T03, T04, and T07) cannot be seen, but the detection can be verified by noting the abundances shown in Table 6. The RGB color composite maps in Fig. 13 show the differences in vegetation modeling for the three methods. The overall greater abundance of grass predicted by the GKLS as compared to FCLS method can be seen by comparing Figs. 13(a) and 13(b). An even greater abundance of grass is predicted by the SSA method, as shown by Fig. 13(c).

The abundance retrieval capability of all methods is good; however, a moderate advantage is shown for the SSA and GKLS methods. Tables 4 and 5 list the model diagnostics for selected points in the scene. The entries in Table 4 list the results for many of the Tyvek panels. The entries in Table 5 list results for selected mixes of vegetation with other backgrounds. Note the RMSE for the GKLS is always less than or equal to the RMSE of the FCLS. This is to be expected since the GKLS algorithm iterates with increasing γ until the RMSE converges to a minimum RMSE. A higher value of γ corresponds to greater nonlinearity in the model. The RMSE values of the SSA method are almost always higher than either of the other methods; however, they were still in an acceptable range. The higher values may be attributed to SSA spectra having overall higher values than reflectance.

Figure 14 shows plots of observed and estimated spectra, as computed by the three methods, for one of the Tyvek panel (TY-5) mixtures with background at location (sample,line)=(X,Y)=(212,76) and one of the vegetation mixtures at (X,Y)=(99,79). Observed spectra (labeled as “O”) are displayed in dark gray and estimated modeled spectra (labeled as “E”) are displayed in lighter gray.

Fig. 14

Plots of observed and estimated spectra are shown for a Tyvek panel mixture with background and a vegetation mixture. Observed spectra (labeled as “O”) are displayed in dark gray and estimated modeled spectra (labeled as “E”) are displayed in lighter gray. Plots for the Tyvek panel T05A mixture with background at location (X,Y)=(212,76) as computed by the FCLS, GKLS, and SSA methods are shown in (a), (b), and (c), respectively. Plots for the vegetation mixture at location (X,Y)=(99,79) are displayed in (d), (e), and (f). Diagnostic Tyvek features are located at 1.72 and 2.31  μm. The shaded areas in the 1.33 to 1.50  μm and 1.79 to 2.01  μm regions indicate where bad bands have been removed because of atmospheric absorption.


Figures 14(a)14(c) show results for the T05A panel mixed with background as computed by the FCLS, GKLS, and SSA methods, respectively. The GKLS result is reported for γ=5.32. Note there are diagnostic Tyvek features located at 1.72 and 2.31  μm (see Fig. 9). Overall, the fits between the observed and estimated spectra are good for the three methods, especially near the diagnostic Tyvek features. Note the GKLS improves the fit as compared to the FCLS method in agreement with the numerical results listed in Table 2.

Figures 14(d)14(f) show results for the vegetation mix as computed by the FCLS, GKLS, and SSA methods. The GKLS result from γ=6.43 is shown. The fits between the observed and estimated spectra were good for the three methods. As with the Tyvek panel, the GKLS improves the fit as compared to the FCLS method in agreement with the numerical results listed in Table 2. However, in this case, the SSA method was slightly better.

Tables 6 and 7 show the abundance results of experiment 2 for the three methods. Table 6 shows these results for the panels. At least one pixel location is listed for each Tyvek panel shown. Although most panels are subpixel, the ground resolution distance of the sensor can cause the panel’s signature to occur over multiple pixels. Not all locations are listed due to an uncertainty factor of the actual sample and also because of space limitations. We notice there is not always a significant difference in abundances between the FCLS and GKLS methods, such as with T04A (one of the pixels in TY-4) and T09B (one of the pixels in TY-9). However, at these locations, the SSA method actually does estimate significantly higher abundances. Overall, the abundances estimated by the SSA method are higher, and for the case of T07 it is sufficiently high to trigger the detection of that panel. The abundance estimates from the other methods applied to T07 are too low to be considered significant.

Table 7 shows the results for vegetation mixtures. The vegetation abundance estimates are notably different for each method. It is difficult to quantitatively assess the results in Table 7 because of a lack of ground truth. However, they corroborate with the qualitative results shown in Fig. 13 and indicate that the SSA method is preferable over the FCLS for better characterizing the nonlinear mixtures of vegetation in the scene.


Results of Experiment 3

The results for experiment 3 (airborne HSI data) are shown in Figs. 15Fig. 1617 and Tables 8 and 9. Figure 15 shows an RGB composite abundance map and five single-layer abundance maps for the fixed-γ GKLS method (γ=5). The composite abundance map visually shows the method’s response to three thicknesses of oil: the R layer depicts the abundance of 4.0-mm-thick oil; the G layer depicts the abundance of 1.85-mm-thick oil; and the B layer depicts the abundance of 1.00-mm-thick oil. The single-layer abundance maps show the results for five of the six endmembers. Due to the small size of these photos in the report, it is difficult to see the detailed patterns, but it is possible to get a sense of the method’s behavior. The oil is expected to spread out in a certain pattern from thick to thin, which is indeed what occurs. In the original photos of the three methods (not shown here but available in supplemental material), there is little difference in the patterns observed between the GKLS and SSA methods. However, there is a notable difference in the patterns seen for the 4.0-mm-thick oil between the FCLS and GKLS (or FCLS and SSA) methods.

Fig. 15

Fraction-plane false-color composites and single-plane abundance maps, using the GKLS method, are shown for five of the six endmembers (water endmember not shown): (a) RGB layer (R=4  mm, G=1.85  mm, R=1.0  mm); (b) single layer (4-mm-thick oil); (c) single layer (1.85-mm-thick oil); (d) single layer (1.0-mm-thick oil); (e) single layer (0.5-mm-thick oil); and (f) single layer (orange oil).


Fig. 16

Abundance maps (fraction planes) of the 4.0-mm-thick oil are shown for the FCLS and GKLS (γ=5) methods: (a) FCLS method: single layer map for 4-mm-thick oil and (b) GKLS method (γ=5): single layer map for 4-mm-thick oil.


Fig. 17

Plots of the observed and estimated spectra at two selected locations: the results for the FCLS, GKLS, and SSA methods at location (x,y)=(401,375) are shown in (a), (b), and (c), respectively. The results for the FCLS, GKLS, and SSA methods at location (x,y)=(323,203) are shown in (d), (e), and (f), respectively. The results for the GKLS methods in (b) and (e) are for γ=5.0. F1 and F2 are diagnostic oil features located at 1.2 and 1.73  μm, respectively. Bad bands caused by atmospheric absorption regions are indicated with the pink-colored strips at 1.10 to 1.16  μm, 1.32 to 1.48  μm, and 1.79 to 1.96  μm.


Table 8.

Abundance results: estimated abundance values of the methods are listed for selected locations (pixels) in the AVIRIS DWH scene.

LocationDWH oilDWH oilDWH oilDWH oil(546, 481)(1154, 112)
Method(x,y)4.0 mm1.85 mm1.0 mm0.5 mmOrangeWater
FCLS(317, 202)0.9650.0000.0000.0000.0000.035
GKLS γ=5(317, 202)0.9560.0000.0000.0000.0000.044
SSA(317, 202)0.9510.0000.0000.0000.0000.049
FCLS(318, 202)0.9880.0000.0000.0000.0000.012
GKLS γ=5(318, 202)0.9730.0000.0000.0000.0000.027
SSA(318, 202)0.9650.0000.0000.0000.0000.035
FCLS(323, 203)0.9660.0430.0000.0000.0090.000
GKLS γ=5(323, 203)0.9560.0440.0000.0000.0000.000
SSA(323, 203)0.9520.0480.0000.0000.0000.000
FCLS(410, 375)0.1330.0000.1710.0000.1870.508
GKLS γ=5(410, 375)0.1370.0000.2300.0520.2010.380
SSA(410, 375)0.1770.0000.0820.3390.1080.294
FCLS(441, 441)0.1190.0000.1420.0000.2930.445
GKLS γ=5(441, 441)0.0730.0000.1810.0940.3030.349
SSA(441, 441)0.0750.0000.0000.4600.1770.288
The bold values indicate the best fitting model for each location.

Table 9

Model diagnostics: RMSE values computed for the fit between the observed and estimated mixed spectra are listed for selected locations of the AVIRIS DWH scene.

(317, 202)0.00600.00640.0274
(318, 202)0.00700.00730.0286
(323, 203)0.00260.00270.0133
(410, 375)0.00920.00730.0287
(441, 441)0.01320.01180.0484

Figure 16 shows the single-layer abundance maps of the 4.0-mm-thick oil for the FCLS and fixed-γ GKLS (γ=5) method. The FCLS method shows much of the same large area of thick oil in the central-right side of the scene as the GKSL (and SSA, not shown). However, these maps also show a notable difference in pattern between the FCLS and GKLS methods, as just mentioned above. The linear FCLS indicates a significant amount of 4.0-mm-thick oil throughout the scene in areas that would not be expected. Physically we expect the oil to spread out in a certain pattern from thick to thin. There is isolated 4.0 mm oil scattered throughout the area in the FCLS photo not seen in the GKLS photo. Neither the GKLS nor the SSA methods indicate much 4.0-mm-thick oil in these regions.

Tables 8 and 9 show the abundance results and RMSE values, respectively, for the three methods at selected locations of the AVIRIS scene. In Table 8, the abundance estimates are not very different for the first three locations. Interestingly, the RMSE values in Table 9 indicate the mixtures are likely to be linear at these locations. The abundance estimates listed in Table 8 for the locations at (x,y)=(410,375) and (x,y)=(441,441) are in agreement with Figs. 15 and 16. The smaller RMSE value for the GKLS at this location indicates the mixture is likely to be nonlinear.

Figure 17 shows plots of the observed and estimated spectra for a thinner layer of oil at (sample,line)=(410,375) and a thicker layer of oil at (sample,line)=(323,203). The fits are qualitatively good except for a significant difference in the F2 diagnostic region for the FCLS method and in the F1 diagnostic region for the SSA method. However, overall the fits are quite good. The larger values of RMSE listed for the SSA method are likely a consequence of the magnitude of the SSA spectra as compared to the normalized reflectance values used in the FCLS and GKLS methods (see range of the y-axis in Fig. 17).



This study has investigated the use of a GKLS method applied to reflectance data, and an SSA method for nonlinear mixture analysis. In the case of the SSA method, an FCLS method is applied to data that have been converted from reflectance space to SSA space. Our baseline method was the FCLS method applied to reflectance data. Our hypothesis is that, for intimate (nonlinear) mixtures, both of these methods will provide improved modeling and abundance estimates as compared to the baseline FCLS method. Three scenarios are tested on different scenes: laboratory, aerial over land, and aerial over oil-contaminated seawater.

Overall the results for experiment 1 using laboratory data from the Pika II sensor indicate the FCLS method has a poor capability for modeling intimate mixtures. In contrast, both the GKLS and SSA methods are much better at retrieving actual material abundances. Whether or not one is better than the other is not conclusive. However, we conclude that our hypothesis is confirmed and that both of these methods provide a better estimate of abundance for mixtures exhibiting nonlinearity. For the laboratory experiment of known abundance quantities, the SSA and GKLS methods responded well to the nonlinearity present in a mixture of materials and provide better estimates of abundance than the linear FCLS method for the DiDy and soda-lime glass bead mixtures.

The GKLS γ parameter determines the degree of nonlinear behavior exhibited by the GKLS method and affects its accuracy for estimating abundances. The automated GKLS method attempts to select the most appropriate gamma based on achieving a minimum of the model’s RMSE. We conclude the RMSE metric seems to respond to a mixture being linear and nonlinear, but is not a reliable metric to determine the degree of nonlinearity present. It could not be used to achieve the most accurate estimate of abundance and an alternative approach should be sought to automate the GKLS method.

For mixtures known to be nonlinear, the fixed-γ implementation of GKLS with γ=5 or 6 provides a good estimate of abundance. The fixed-γ GKLS and the SSA methods can be computed in approximately the same amount of time and provide approximately the same accuracy for estimating abundances. Compared to the fixed-γ GKLS method, the automated GKLS was much slower to compute by a factor of about 19 and did not really achieve better accuracy.

For experiments 2 and 3 using the airborne data from HyMap (scene over land) and AVIRIS (scene over water), we also conclude that our hypothesis is confirmed and that either of these methods provides a better estimate of abundance for mixtures exhibiting nonlinearity. Unfortunately, the data set used in experiment 2 was from a collection campaign not particularly focused on nonlinear mixing problems, and although it contained many nonlinear mixtures, precise ground truth is not available. However, the data provide good opportunity to observe that linear mixtures did not significantly degrade the performance of the nonlinear GKLS and SSA methods. In experiment 3 (Deepwater Horizon oil spill), the SSA and GKLS methods performed equivalently and appeared to respond better to the nonlinear mixtures of oil thicknesses in the water than the linear method. Specifically, the linear method appears to be predicting too much of the thicker oils. This can only be assessed qualitatively because, here, too, there is a lack of ground truth. Additional analysis should be performed on scenes containing a variety of known nonlinear mixtures that are well documented with ground truth.

Appendix provides information about obtaining data from the first experiment. All of the data as well as dark field data cubes (with 80, 120, and 240 bands) and some additional information [e.g., the ENVI ROI files and an ENVI ‘Spectral Math’ expression file for implementing Eqs. (2) and (3)] are available for download as one compressed (.zip) file.

Additional experiments are recommended that investigate the effects of adding endmembers to the models and also to the physical mixtures in a controlled environment. In experiment 1, the number of controlled endmembers was limited to two.



Obtaining the Data

The focus of the first experiment discussed in this study has been on the analysis of a single Pika II HSI data cube (file: 20130521_sphere_array_2.bil). However, 15 measurements, in total, of the glass beads in the 96-well sample plate were acquired as listed in Table 10. Two are simply redundant measurements (640  samples×500  lines×80  bands prior to spectral subsetting in ENVI). Six others are of the same sample plate and glass bead mixtures but with 120 and 240 bands. These data are somewhat noisier than those with 80 bands but were acquired to assess the impact of an increased noise level on abundance estimations. Three cubes are with the sensor moved closer to the sample plate (still with the Xenoplan objective lens) and with 80 bands. Two additional cubes contain two replicate sets of glass bead mixtures placed in adjacent empty cells of the plate (for a total of 15 filled cells). These cubes, too, contain 80 bands. All of the data sets are with four QTH lamps and nadir-looking sensor configuration shown in Fig. 1. To date, no data have yet been acquired with a line-light source to yield true bidirectional reflectance spectra. All of these data sets as well as dark field data cubes (with 80, 120, and 240 bands) and some additional information [e.g., the ENVI ROI files and an ENVI “spectral math” expression file for implementing Eqs. (2) and (3)] are available for download as one compressed (.zip) file. Send an e-mail to requesting the data; a link for downloading will be sent in return. All questions and comments about these data should be addressed to R.G. Resmini.

Table 10

A list is shown of the available Pika II data sets of the glass beads. (File 20130521_sphere_array_shifted_8.bil is the same as file 20130521_sphere_array_8.bil except it has been shifted so that some regions of interest defined in other cubes in the data set also apply.)



References are made to certain commercially available products in this paper to adequately specify the experimental procedures involved. Such identification does not imply a recommendation or endorsement by the U.S. government, nor does it imply that these products are the best for the purpose specified. Thanks are given to Dr. Charles Bachmann and the U.S. Naval Research Laboratory for providing the HyMap imagery, as well as associated aerial photographs and atmospheric calibration for experiment 2. Thanks are also given to Dr. Roger Clark and Eric Livo (USGS) for the data provided in experiment 3, including the laboratory and field measurements, as well as the AVIRIS imagery. The MITRE Innovation Program (MIP) is gratefully acknowledged for funding the HSI microscopy aspect of the project in which the study presented here was conducted. This article has been approved by NGA for public release, Case 16-132.



J. Adams, M. Smith and P. Johnson, “Spectral mixture modeling: a new analysis of rock and soil types at the Viking Lander 1 site,” J. Geophys. Res., 91 8098 –8112 (1986). Google Scholar


J. Boardman, “Automating linear mixture analysis of imaging spectrometry data,” in Proc. of the Int. Symp. on Spectral Sensing Research (ISSSR), (1994). Google Scholar


R. S. Rand, “A physically-constrained localized linear mixing model for TERCAT applications,” Proc. SPIE, 5093 547 (2003). PSISDG 0277-786X Google Scholar


R. S. Rand, “Automated classification of built-up areas using neural networks and subpixel demixing methods on multispectral/hyperspectral data,” in Proc. of the 23rd Annual Conf. of the Remote Sensing Society (RSS97), (1997). Google Scholar


R. S. Rand, “Exploitation of hyperspectral data using discriminants and constrained linear subpixel demixing to perform automated material identification,” in Proc. of the Int. Symp. on Spectral Sensing Research (ISSSR), (1995). Google Scholar


R. S. Rand and D. M. Keenan, “A spectral mixture process conditioned by Gibbs-based partitioning,” IEEE Trans. Geosci. Remote Sens., 39 (7), 1421 –1434 (2001). Google Scholar


D. C. Heinz and C.-I. Chang, “Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., 39 (3), 529 –545 (2001). IGRSD2 0196-2892 Google Scholar


D. Montgomery and E. Peck, “Introduction to linear regression analysis,” Wiley Series in Probability and Mathematical Statistics, 2nd ed.John Wiley and Sons, Hoboken, New Jersey (1992). Google Scholar


M. Eismann and D. Stein, “Stochastic mixture modeling,” Hyperspectral Data Exploitation: Theory and Applications, Wiley, Hoboken, New Jersey (2007). Google Scholar


O. Eches et al., “Bayesian estimation of linear mixtures using the normal compositional model–application to hyperspectral imagery,” IEEE Trans. Image Process., 19 (6), 1403 –1413 (2010). Google Scholar


L. Zhuang et al., “Normal endmember spectral unmixing methods for hyperspectral imagery,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 8 (6), 2598 –2606 (2015). Google Scholar


A. Castrodad et al., “Sparse modeling for hyperspectral imagery with lidar data fusion for subpixel mapping,” in IEEE Int. Geoscience and Remote Sensing Symp., (2012). Google Scholar


A. Castrodad et al., “Learning discriminant sparse representations for modeling, source separation, and mapping of hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., 49 (11), 4263 –4281 (2011). IGRSD2 0196-2892 Google Scholar


J. Bhatt, J. Manjunath and M. Raval, “A data driven stochastic approach for unmixing hyperspectral imagery,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 7 (6)), 1936 –1946 (2014). Google Scholar


Y. Yuan, Y. Feng and X. Lu, “Projection-based NMF for hyperspectral unmixing,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 8 (6), 2632 –2643 (2015). Google Scholar


W. Wang and Y. Qian, “Adaptive L1/2 sparsity-constrained NMF with half-thresholding algorithm for hyperspectral unmixing,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 8 (6), 2618 –2631 (2015). Google Scholar


N. Keshava and J. F. Mustard, “Spectral unmixing,” IEEE Signal Process. Mag., 19 (1), 44 –57 (2002). ISPRE6 1053-5888 Google Scholar


B. Hapke, Theory of Reflectance and Emittance Spectroscopy, 455 Cambridge University Press, New York (1993). Google Scholar


J. F. Mustard and C. M. Pieters, “Photometric phase functions of common geologic minerals and application to quantitative analysis of mineral mixture reflectance spectra,” J. Geophys. Res., 94 (B10), 13619 –13634 (1989). JGREA2 0148-0227 Google Scholar


J. F. Mustard and C. M. Pieters, “Quantitative abundance estimates from bidirectional reflectance measurements,” J. Geophys. Res., 92 (B4), E617 –E626 (1987). JGREA2 0148-0227 Google Scholar


J. F. Mustard and C. M. Pieters, “Abundance and distribution of ultramafic microbreccia in moses rock dike: quantitative application of mapping spectroscopy,” J. Geophys. Res., 92 (B10), 10376 –10390 (1987). JGREA2 0148-0227 Google Scholar


S. G. Herzog and J. F. Mustard, “Reflectance spectra of five component mineral mixtures: implications for mixture modeling,” Lunar Planet. Sci., 27 535 –536 (1996). Google Scholar


R. G. Resmini et al., “Constrained energy minimization applied to apparent reflectance and single-scattering albedo spectra: a comparison,” Proc. SPIE, 2821 3 –13 (1996). Google Scholar


R. G. Resmini, “Enhanced detection of objects in shade using a single-scattering albedo transformation applied to airborne imaging spectrometer data,” in The Int. Symp. on Spectral Sensing Research, 7 (1997). Google Scholar


J. M. P. Nascimento and J. M. Bioucas-Dias, “Unmixing hyperspectral intimate mixtures,” Proc. SPIE, 7830 8 (2010). PSISDG 0277-786X Google Scholar


N. Dobigeon et al., “Nonlinear unmixing of hyperspectral images: models and algorithms,” IEEE Signal Process. Mag., 31 (1), 82 –94 (2014). ISPRE6 1053-5888 Google Scholar


R. Close et al., “Using physics-based macroscopic and microscopic mixture models for hyperspectral pixel unmixing,” Proc. SPIE, 8390 83901L (2012). Google Scholar


R. Heylen, M. Parente and P. Gader, “A review of nonlinear hyperspectral unmixing methods,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 7 (6), 1844 –1868 (2014). Google Scholar


J. Li, X. Li and L. Zhao, “Nonlinear hyperspectral unmixing based on sparse non-negative matrix factorization,” J. Appl. Remote Sens., 10 (1), 015003 (2016). Google Scholar


A. Marinoni et al., “Nonlinear hyperspectral unmixing using nonlinearity order estimation and polytope decomposition,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 8 (6), 2644 –2654 (2015). Google Scholar


A. Marinoni, A. Plaza and P. Gamba, “Harmonic mixture modeling for efficient nonlinear hyperspectral unmixing,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., 9 (9), 4247 –4256 (2016). Google Scholar


H. Kwon and N. M. Nasrabadi, “Kernel matched subspace detectors for hyperspectral target detection,” IEEE Trans. Pattern Anal. Mach. Intell., 28 (2), 178 –194 (2006). ITPIDJ 0162-8828 Google Scholar


G. Camps-Valls and L. Bruzzone, “Kernel-based methods for hyperspectral image classification,” IEEE Trans. Geosci. Remote Sens., 43 (6), 1351 –1362 (2005). IGRSD2 0196-2892 Google Scholar


B. Scholkopf and A. J. Smola, Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond, The MIT Press, Cambridge, Massachusetts (2002). Google Scholar


J. B. Broadwater et al., “Kernel fully constrained least squares abundance estimates,” in Proc. of the IEEE Int. Geoscience and Remote Symp. (IGARSS 2007), 4091 –4044 (2007). Google Scholar


J. B. Broadwater and A. Banerjee, “A comparison of kernel functions for intimate mixture models,” in Proc. of the IEEE First Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), (2009). Google Scholar


J. B. Broadwater and A. Banerjee, “A generalized kernel for areal and intimate mixtures,” in Proc. of the IEEE First Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), (2010). Google Scholar


J. B. Broadwater and A. Banerjee, “Mapping intimate mixtures using an adaptive kernel-based technique,” in Proc. of the IEEE First Workshop on Hyperspectral Image and Signal Processing: Evolution in Remote Sensing (WHISPERS), (2011). Google Scholar


R. S. Rand, A. Banerjee and J. Broadwater, “Automated endmember determination and adaptive spectral mixture analysis using kernel methods,” Proc. SPIE, 8870 88700Q (2013). PSISDG 0277-786X Google Scholar


R. G. Resmini et al., “An analysis of the nonlinear spectral mixing of didymium and soda lime glass beads using hyperspectral imagery (HSI) microscopy,” Proc. SPIE, 9088 9088OZ (2014). PSISDG 0277-786X Google Scholar


R. N. Clark et al., “A Method for quantitative mapping of thick oil spills using imaging spectroscopy,” (2010). Google Scholar


R. S. Rand, R. N. Clark and K. E. Livo, “Feature-based and statistical methods for analyzing the Deepwater Horizon oil spill with AVIRIS imagery,” Proc. SPIE, 8158 81580N (2012). Google Scholar


J. R. Schott, Remote Sensing: The Image Chain Approach, 688 2nd ed.Oxford University Press, New York (2007). Google Scholar


R. P. Brent, Algorithms for Minimization without Derivatives, Prentice-Hall, Englewood Cliffs, New Jersey (1973). Google Scholar


, “Pika II – VNIR hyperspectral imaging camera,” (2016) December 2013). Google Scholar


B.-C. Gao et al., “Atmospheric correction algorithm for hyperspectral remote sensing of ocean color from space,” Appl. Opt., 39 (6), 887 (2000). Google Scholar


B.-C. Gao and C. O. Davis, “Development of a line-by-line-based atmosphere removal algorithm for airborne and spaceborne imaging spectrometers,” Proc. SPIE, 3118 132 –141 (1997). Google Scholar


Robert S. Rand received his BS degree in physics from the University of Massachusetts at Lowell and his PhD in engineering physics from the University of Virginia, Charlottesville. He has worked for agencies with the U.S. Department of Defense since 1977. His recent efforts have focused on hyperspectral exploitation using feature-based methods and spectral mixture models, as well as the use of three-dimensional information from LIDAR to improve the exploitation of hyperspectral imagery.

Ronald G. Resmini received his MS degree in geology from Boston College, Boston, Massachusetts, USA, and his PhD in geology from Johns Hopkins University, Baltimore, Maryland, USA. He is a research scientist at the MITRE Corporation, McLean, Virginia, USA, and an associate professor in the College of Science at George Mason University, Fairfax, Virginia, USA. He specializes in visible to infrared multi- and hyperspectral remote sensing, the geological and geophysical sciences, and algorithm development for remote sensing applications.

David W. Allen received his PhD in earth systems and geoinformation sciences with a concentration in remote sensing from George Mason University. He is a research chemist at the NIST, where he has worked for over 20 years. Current research interests include advancing hyperspectral imaging for applications related to medicine, security, and the environment. This effort encompasses the development of end-to-end analysis methodology that integrates standards and best practices for hyperspectral imager performance metrics.

© The Authors. Published by SPIE under a Creative Commons Attribution 3.0 Unported License. Distribution or reproduction of this work in whole or in part requires full attribution of the original publication, including its DOI.
Robert S. Rand, Ronald G. Resmini, and David W. Allen "Modeling linear and intimate mixtures of materials in hyperspectral imagery with single-scattering albedo and kernel approaches," Journal of Applied Remote Sensing 11(1), 016005 (11 January 2017).
Received: 19 July 2016; Accepted: 13 December 2016; Published: 11 January 2017

Back to Top