2.1.

Data

The baseline image is a 49-band LWIR hyperspectral radiance image spanning 7.93 to $12.56\mu \mathrm{m}$ window. The material display board was imaged under a cold sky. The exact geometry of the setup is not known; the board was angled at 45 deg relative to the surface normal. The viewing geometry was slightly downward looking at an approximate 10-deg angle from a high rise building at an unknown distance. The sensor had a slight focusing issue during the data collection. Readers can see a spectral radiance plot of calcite taken from this dataset in Fig. 2 along with a three-channel image of the material display board.

Fig. 2

(a) A calcite emissivity spectrum is plotted with data from a scene pixel (in radiance) that contains calcite (see red box in Fig. (b)]. (b) An RGB image was created from three hyperspectral bands to give the reader context of the material display board. The squares are blocks of various materials. The warm and emissive materials show up as bright white squares; metallic materials are reflecting cold sky and appear as black squares; and the bottom row shows several purple squares, these materials have spectral features located over the three hyperspectral bands used to make this image.

The measured radiance is affected by the material spectral emissivity, bidirectional reflectance distribution, geometry, and temperature. Atmospheric radiance and transmission also affect the measured radiance; these factors are influenced by the atmospheric gas’ transmission, pressure, number density, concentration, and temperature.

In this work, an assumption is made that all materials are diffuse reflectors. The Lambertian (diffuse) LWIR radiance model² is defined as

The target spectra came from the Arizona State University Spectral Library and the Jet Propulsion Lab EcoStress Spectral Library (the original measurement was made by USGS). A plot of the target spectra is shown in Fig. 3:

Fig. 3

The target material spectra used in this study have a diverse set of spectral features of varying widths and depths.

2.2.

Band Selection and Multispectral Band Creation

Groups of contiguous bands are selected from the hyperspectral radiance cube. The groups are selected such that the edge of the adjacent group does not overlap. Each group is then averaged spectrally. The resulting set is intended to simulate data collected by a multispectral system. Figure 4 provides an example of a four-band case.

Fig. 4

The original dataset is a full hyperspectral cube, groups of contiguous bands are then selected from this cube, and these groups are averaged to create the individual multispectral bands. “$N_x$” defines the number of HSI bands in the group.

It should be noted that detector noise characteristics for the multispectral system were not simulated.

2.3.

Data Normalization by Brightness Temperature

To use the multispectral data in a target detection algorithm, a normalization technique was used as an initial processing step to help reduce variability caused by the material temperature. This is done by first computing the brightness temperature ($T_{\max }$ [K]) of every multispectral pixel radiance vector $im(p)$:

The brightness temperature is then used to normalize the data by dividing each pixel by a Planck curve defined by the brightness temperature:

2.4.

ACE

The ACE algorithm is a popular hyperspectral target detection algorithm.^2,6,12 The algorithm uses the scene data multivariate mean and covariance in a process referred to as whitening to suppress the data’s direction of maximum variability. The ACE algorithm whitens and normalizes the pixel vectors ($x$) and target vector(s) ($t$) and then computes cosine angle between the two. The target spectra are defined in Sec. 2.1. ACE is defined as

2.5.

Separation Score

Ideally, the target detection algorithm will have a good rate of true positive detection and a low rate of false positive detection. While it is not exactly analogous, if the algorithm is performing well, the output scores over target areas will be higher than those over background areas. The assumption used here is that a greater separation in ACE scores between target and background will result in an improvement in detection performance (either a decrease in false positive rate or an increase in true positive rate).

As mentioned in Sec. 1.1, the sensor experienced a spatial focusing issue which resulted in some pixels being mixed with both target and background. In this work, the number of on-target pixels was much less than the number of background pixels. Figure 5(a) shows a detection plane for the calcite target. Plotting a histogram of the ACE scores for on-target and background regions, we see these two groups overlap where the pixels are mixed with either target or background [see Fig. 5(b)].

Fig. 5

(a) Observe the ACE detection plane for the calcite target. The high detection values (bright yellow square) span the target region. The blue/green region is entirely background. (b) The histogram of ACE scores for the target area and background area shows perfect separation between the pure target and background areas. Pixels that are mixed with target and background overlap either class. The $y$-axis uses a log scale to account for the significant imbalance between the pixel counts of the two classes.

There is a relatively small number of target pixels and many of the ACE scores in the target-region are statistical outliers. For this reason, a median is used to characterize the target-region ACE scores. The background region, which contains a much greater number of pixels, was characterized by taking a mean value for this region. Capturing the amount of separability between the target and background is done by differencing the median of the target area ACE scores and the mean of the background ACE scores:

2.6.

Parameterization

Parameterizing this problem is done by defining the lower and upper bounds of each multispectral bandpass. An illustration of this is provided for the four-band case in Fig. 6. The multispectral system’s bandpasses can be defined by a vector $Z=(z_1,z_2,\dots ,z_N)$, which is a vector of hyperspectral band indices. N.B. the $N$ defined here is $2\times$ the number of multispectral bands.

Fig. 6

An illustration of the four-band case. The blue rectangles represent the bandpass filters. The bandpass filters are constrained to have some width, potentially contiguous but not overlap.

There are three constraints that are used to define this problem. Constraint 1: Vector $Z\in {\mathbb{Z}}^{+}$, that is, it is an integer vector since it contains hyperspectral band indices. Constraint 2: Vector $Z_n\in [1,49]$, the HSI used here has 49 bands and therefore the band indices ($n$) must be between 1 and 49. Constraint 3: $z_1<z_2,z_2\le z_3,\dots z_{N-1}<z_N$, the multispectral bandpasses must have a width of at least one band and they are allowed to be contiguous. Parameterizing Eq. (5), we can define the separation $\mathrm{Score}$ for a given scenario, $Z$, as:

4.1.

Integer Particle Swarm Optimization

Particle swarm optimization (PSO) interrogates a problem with a population (or swarm) of individual potential solutions (particles). As the algorithm makes iterative steps toward an optimal solution, it is influenced by the previous global best solution as well the individual particle’s best previous solution. Each particle has a velocity that influences its momentum as it traverses the solution space. Having momentum as well as influence from the global best solution helps prevent particles from getting stuck in local minima.

The implementation used in this work closely follows¹⁴ which is similar to continuous PSO; however, the dependent variables are rounded to their nearest integer and the constraints mentioned in Sec. 2.6 are enforced. The optimization was setup to use: 50 particles, 200 iterations ($M=200$), a cognitive weight $(c_1)$ of 0.6, a social weight $(c_2)$ of 0.6, a velocity constriction factor $(\chi )$ of 0.729, inertia weight ($\omega$) of 0.9, and an inertial weight decay function of:

4.2.

Dual Annealing

Simulated annealing is a well-known heuristic for unconstrained optimization problems. It is an improvement over basic local search methods because it allows the incumbent solution to move in a direction that is not locally improving with a certain probability, thus allowing the solution to escape local optima. Dual annealing is a type of simulated annealing algorithm that combines classical simulated annealing and fast simulated annealing. Dual annealing is a standard routine in the SciPy Optimize package.^15,16

Simulated annealing routines are generally applied to unconstrained optimization problems. Since the problem has some constraints—all bands have lower and upper bounds, bands cannot overlap—the objective function was modified to return a poor score if one of the constraints was violated. In the experiments, the total number of iterations was limited to 300, and the visit parameter was set to 2.9 (instead of the default value of 2.62). This allows the algorithm to jump to greater distances in the search space. The remaining parameters were set to their default values in SciPy Optimize.

4.3.

Differential Evolution

Differential evolution is an evolutionary algorithm. Starting with a randomly selected set of candidate solutions, the algorithm iterates through a sequence of mutation, recombination, and selection steps until a termination criterion is met. The algorithm is also included in the SciPy Optimize package.¹⁷

The SciPy Optimize package implementation of differential evolution allows constraints to be specified. We added constraints to limit the upper and lower bounds of the bands and constraints to ensure that the bands do not overlap. We set the maximum number of iterations to 300 and used default values in SciPy Optimize for all other parameters.

5.1.

Brute Force Processing

The brute force global optimum separation scores and associated bandpasses for the three- and four-band scenarios are provided in Table 1 of the Appendix. The bandpasses are also displayed graphically in Figs. 7 and 8. While the separation score is a useful metric for the optimization algorithm, detection performance is best captured by examining the probability of detection (PD) and the probability of false alarm (PFA). This done using a receiver operator characteristic (ROC) curve. Overall performance can be quantified by taking the area under the ROC (AUROC) curve. Figures 7 and 8 show ROC curves for the three- and four-band cases. Figure 7 shows expected behavior for calcite and gypsum, where the optimal spectral band selection is aligned with dominant spectral features. However, this does not occur with the individual case of limestone or for the three-material case used in the optimization. It is unclear why this happened. In the individual limestone case, the other subtle spectral features (8.5, 9.2, and $9.8\mu \mathrm{m}$) may have offered more information relative to the background spectra than the dominate feature ($11\mu \mathrm{m}$). In the three-material case, one material may have gained substantially from the selected band configuration and mask performance gains (or losses) for the other materials. Figure 8 shows more consistent behavior for the four-band case, where bands fall on dominant spectral features. Examining the ROC curve, in both cases, the contiguous multispectral bands perform worse than the hyperspectral system. An exciting result is that a multispectral system with band optimization for three materials can produce better performance than the original hyperspectral system. Optimizing for each individual material results in further improvement for the material under inspection; however, detection performance for the other materials declines substantially. It should be noted that this observation applies to only this dataset. Another dataset with the same target materials, but different background, could result in a different optimal band combination. Future work could examine the extensibility of these bandpasses, comparing these results to results collected over different backgrounds or different times of day (different thermal conditions).

Fig. 7

The brute force result for the three-band case is shown here. Each vertical pair of plots corresponds to an optimization scenario. For each plot pair, the upper plot shows the three material spectra of interest and selected spectral bandpasses. The lower plot shows the ROC curves and AUROC curve values for each of the three materials of interest. This shows a comparison between the original hyperspectral data, three-band contiguous data, and three-band optimized data.

Fig. 8

The brute force results from the four-band case.

5.2.

Numerical Optimization Results

The detailed optimization results are provided in Table 1 of the Appendix. For each scenario, 30 optimization runs were performed. The mean, standard deviation, median, minimum, and maximum values are given for these 30 optimizations. The bandpasses for the maximum separation score are also given.

Figure 9 provides a summary of algorithm performance. The mean is calculated for the four material types of a single band number and optimization type. For example, the mean of the three-band particle swarm algorithm 30-run means is mean (0.518, 0.576, 0.886, 0.504) = 0.621. And the mean of the 30-run standard deviations is mean (0.031, 0.105, 0.019, 0.028) = 0.046.

Fig. 9

Bar chart shows a summary of the separation scores for each algorithm. The bar values are produced by taking the mean value of the mean separation scores for the four material types for each algorithm. The error bars are produced by taking the mean value of the standard deviation of the separation scores. The dual annealing and differential evolution algorithms produced the highest and most consistent results for the three-band case. For the four- and five-band cases, all algorithms are similar. It is worth noting that the particle swarm algorithm has higher variance than the other two algorithms, which is a good thing if the optimization can be run many times—as a higher maximum score could be more attainable with this algorithm.

For the three-band case, all three algorithms performed well and found the global maximum for several materials. The particle swarm algorithm found the global maximum for all materials; however, the separation score standard deviation was highest for all materials using this algorithm making it the most inconsistent performer. The differential evolution algorithm produced the best results on average with a higher mean separation score and a lower standard deviation.

In the four-band case, each algorithm found the global maximum in two out of the four material scenarios. The differential evolution algorithm had the highest mean separation score and the lowest mean standard deviation for the 30-run means.

The five-band case has many more potential combinations, and it was therefore not possible to compute a brute force result. The particle swarm algorithm produced the best results on average; however, there was large variability in the scores. The dual annealing algorithm produced the most consistent results but only by a small margin compared with the differential evolution. The performance margins here are closer than in the three-band case.

Figure 10 shows the best five-band optimization results for each material. The bandpasses for the individual materials follow the dominant spectral features. However, the three-material optimization lacks a bandpass over the dominant limestone feature. If the bandpasses are modified to include a band over the dominant limestone feature an interesting observation is made (see “User Defined” plot in Fig. 10). The AUROC values for limestone increase while the values of calcite and gypsum both decrease. On average, the AUROC values for the optimized bands are higher than the user defined bands. Three-material AUROC average: (0.9879 + 0.9907 + 0.9484)/3 = 0.9757, user defined AUROC average: (0.9569 + 0.9854 + 0.96168)/3 = 0.9680. The optimization algorithm would have likely favored improved performance for calcite and gypsum at the expense of the limestone performance.

Fig. 10

Five-band results obtained from the best numerical optimizations, along with a “User Defined” set of bandpasses in the upper left.

6.1.

Brute Force Computing

A brute force search of all possible spectral band combinations was performed for a three- and four-band system using a 49-band HSI as a data source. Executing this task required using 128 nodes on LLNL’s Syrah supercomputer and 1 week to complete. Given the successful performance of the optimization algorithms, there is little to be gained by brute force searches of this kind in future work; the resources required are significantly greater than what is used for the numerical optimization approach.

6.2.

Numerical Optimization

Integer programing, which is a type of numerical optimization that utilizes integer valued variables, was used to select band combinations from the 49-band HSI. Three algorithms were used in the optimization: particle swarm, dual annealing, and differential evolution. Overall, the performance for all three algorithms is similar. There are many instances where these algorithms found the global optimum solution. The differential evolution algorithm produced high separation scores more consistently than the other algorithms but only by a small margin. The particle swarm algorithm produced a higher standard deviation of separation scores, making it the most inconsistent performer. That said, it did achieve some of the best scores; most notably for the three-band case, where the global maximum was found for all materials.

Plotting the optimized bandpasses against the target vectors revealed that the results were not always intuitive; there are some instances of dominant spectral features not being selected for bandpasses. The most unusual example of this being the three-band limestone optimization where smaller features were selected over the dominant feature. There are several factors that could have caused this, particularly the scene whitening step that is used in ACE which can be heavily influenced by the background. This behavior was also observed in the three-material optimizations for three bands and five bands. When three materials are used in the optimization, it is possible that substantial performance gains of one material might overshadow the decreased performance of the other two.

6.3.

User Defined Optimization

In Fig. 10, the authors modified the three-material optimization by moving one band to the dominant (and unused) limestone spectral feature. The band that was moved existed beyond $12\mu \mathrm{m}$ in a location that was spectrally flat for the targets used here. Moving this band improved performance for limestone detection but resulted in decreased performance for calcite and gypsum. This serves as a reminder that flat target spectra can still be informative relative to the background spectral shape. In this case, the best detection performance overall was obtained by ignoring the primary limestone feature.

6.4.

Future Work

Given the good optimization performance, these tests should be expanded to datasets collected under different thermal conditions, different backgrounds, and single pixel detection. Future work might also include modifying the objective function to penalize optimization steps that result in large performance gains for one material at the expense of other materials.