A Semiparametric Model for Hyperspectral Anomaly Detection

Abstract

Using hyperspectral (HS) technology, this paper introduces an autonomous scene anomaly detection approach based on the asymptotic behavior of a semiparametric model under a multisample testing and minimum-order statistic scheme. Scene anomaly detection has a wide range of use in remote sensing applications, requiring no specific material signatures. Uniqueness of the approach includes the following: (i) only a small fraction of the HS cube is required to characterize the unknown clutter background, while existing global anomaly detectors require the entire cube; (ii) the utility of a semiparematric model, where underlying distributions of spectra are not assumed to be known but related through an exponential function; (iii) derivation of the asymptotic cumulative probability of the approach making mistakes, allowing the user some control of probabilistic errors. Results using real HS data are promising for autonomous manmade object detection in difficult natural clutter backgrounds from two viewing perspectives: nadir and forward looking.

1. Introduction

Hyperspectral (HS) sensors collect the radiation over a wide range of contiguous spectral bands, with each band corresponding to a unique spectral value. The field of view of the sensor is broken into hundreds of thousands of pixels, with each pixel representing from less than one to many squared meters of the region of interest depending on the spatial resolution of the sensor and the height of the sensor during the data collection. A collection of spatial-spectral images is put together resulting in an HS data cube, where the length and width represent the spatial dimension, and the depth represents the spectral dimension [1]. The resulting HS data cube consists of hundreds of thousands of pixels. Each pixel has tens or hundreds of data points, each point corresponding to a unique spectral value. In theory, the spectral signature of each pixel should uniquely characterize the physical material in that spatial land area. In practice, the recorded spectral signatures will never be identical for samples of the same material. Owing to the different illumination conditions, atmospheric effects, sensor noise, and so forth, the resulting spectral signatures for HS imagery pixels containing similar materials will exhibit spectral variability.

The discriminant capability, however, of spectral signatures has led to two major applications: object classification and target detection. The former aims to assign all pixels of the image to thematic classes, the latter searches the image for the presence of specific material (the target). As highlighted in [2], from a theoretical point of view, target detection can be considered as a binary classification problem, aiming at classifying the image into the target class and the background class. But, since targets are usually sparsely populated in the scene, while the background is abundant and heterogeneous, a major distinction between the approaches designed for target detection and classification is that target detectors cannot use criteria based on the minimization of classification error since that would lead to labeling every pixel as background. So, a typical solution proposed in the literature for target detection is to use the Neyman-Pearson approach, as discussed by Manolakis in [3], where maximizing the detection probability for a fixed false-alarm rate is the adopted criterion for the algorithm design.

Due to the availability of spectral signature libraries for a wide range of materials, spectral signature-based target detectors have been widely examined [3, 4]. These methods assume the target spectral signature is both reliable and known a priori and aim at finding highly correlated spectra in the scene corresponding to the reference target spectra; these methods are also known as target matching.

Target matching approaches are complicated by the large number of possible objects of interest and uncertainty as to the reflectance/emission spectra of these objects. For example, the surface of an object of interest may consist of several materials, and the spectra may be affected by weathering or other unknown factors. One may be interested in a large number of possible objects each with several signatures. Thus, the multiplicity of possible spectra associated with the objects of interest and the complications of atmospheric compensation, as well spectral calibration, acquisition geometry, and contamination from adjacent objects (see, for instance, the discussion in [5, 6]), have led to the development and application of anomaly detectors that seek to distinguish observations of unusual materials from typical background materials without reference to target signatures.

Anomalies are defined with reference to a model of the background. Typically, background models are developed adaptively using reference data from either a local neighborhood of the test pixel or a large section of the image. Local and global spectral anomalies are defined as observations that deviate in some way from the neighboring clutter background or the image-wide clutter background, respectively. Both approaches have their merits.

Local anomaly detectors are typically designed under the assumption that an anomalous material (the target) is spectrally distinct from local neighborhood spectra, which are assumed to be controlled by a known multivariate distribution (Gaussian); also, it is assumed that the scales of targets are known a priori, or the viewing perspective is assumed to be nadir and the altitude of the flying platform carrying the sensor is available for target scale estimation. This kind of detectors is susceptible to unknown spectral mixtures that are often obtained by sampling spectra through a moving window in the imagery, as the window is placed at a transition of distinct regions, forcing the neighboring spectral mixture to be compared against spectra of one of the regions in that transition; this may significantly increase the false alarm rate. The local anomaly detector is also susceptible to false alarms that are isolated spectral anomalies. For example, consider a scene containing isolated trees on a grass plain. Each separate tree may be detected as a local spectral anomaly even if the image contains a separate region with many pixels of trees. The most popular local anomaly detector in the HS research community is based on maximum likelihood estimation under the multivariate normal distribution; this detector is commonly known as the Reed-Xi (RX) algorithm [7] and has become a benchmark for comparison. Kwon and Nasrabadi proposed a kernelized version of RX in [8], and Matteoli et al. proposed a segmented version of RX in [9], for the estimation of the local background covariance matrix from global background statistics to cope with the small sample size problem in estimating covariance from local patches in the image; a small sample size problem occurs when the number of spectral observations is lesser than the number of spectral bands (see examination of this problem in [2]). Banerjee et al. in [10] leveraged the employment of kernel methods and the method known as support vector data description (SVDD) to propose the kernel SVDD. But, it is widely known that the performance of kernel methods crucially depends on the kernel function and its parameter(s) [11]. More recently, Gurram and Kwon in [12] and Khazai et al. in [13] have also addressed the parameter sensitivity of kernel-SVDD based detectors, which is still an open area of research. In the statistical based arena, Stein et al. in [14] presented an overview of the statistical anomaly detectors derived from three background models: local Gaussian model, Gaussian mixture model (GMM), and linear-mixing model. More recently, these models were compared against others approaches using the same dataset [15, 16]. In [15], the algorithms RX, GMM, and a cluster-based one were examined. Matteoli et al. in [16] extended the comparison to include the orthogonal subspace projection (OSP) detector and a deterministic signal subspace processing detector. Other classic approaches have also been adapted to the local anomaly detection problem, for example, Fisher’s linear discriminant (see an implementation in [17]).

Global anomaly detectors are based on a simple universal distribution, which is designed to represent the background process in the whole image, thus, the name global. Example of these detectors is the GMM [14] or a different version of the RX algorithm, which estimates its required parameters (mean and covariance) not locally, but using spectra from the entire HS data cube; this version is informally known as the Global RX. By design, these methods are more resistant to the small sample size problem mentioned earlier, but they have limited ability to adapt to all nuances of the background class (sometimes referred to as an underfitting problem), which may result in both high false alarm and low anomaly detection rates. Other earlier versions of global detectors require that an HS data cube is first segmented into its constituent material classes, so detection is achieved by applying a cutoff threshold and automatically locating pixel clusters with pixel values above the threshold, representing the outliers of these classes. These hybrid algorithms vary in the method of segmentation, but also tend to use maximum likelihood detection under the multivariate normal distribution. The stochastic expectation maximization clustering algorithm by Stocker in [18] is a related example; see also Masson and Pieczynski in [19]. But, since segmentation results are known to be also sensitive to algorithms’ parameters (see, for instance, [1]), utility of segmentation algorithms in the context of anomaly detection has not met expectations for varying real world scenarios.

Independently of the kind of anomaly detector in use, the following is a key consideration that should not be ignored: susceptibility to unknown spectral mixtures of unknown distributions often observed by sampling spectra through a moving window in the imagery, where the spectra in test belong to one of the components in that mixture (for instance, a local patch of canopy being compared against neighboring spectral mixture of the same canopy type and a patch of soil). Rosario in [20] examined this particular problem using near infrared HS imagery, where he showed that even on simple real case scenarios (e.g., motor vehicles parked at an open grassy field having also trees in the background on a sunny day) transitions of distinct regions may contribute to 18% of all of the locally observed spectra in the imagery, using a small moving window (10 by 10 pixels). Of course, with increased scene complexity (increased heterogeneity), this percentage reaches higher levels.

This paper presents a quasiglobal, semiparametric anomaly detection (QG-SemiP) approach that is less susceptible to problems and issues mentioned above for existing local or global anomaly detectors. In particular, the approach requires only a small fraction of the HS cube to characterize the unknown clutter background (hence, the term quasi-global), in contrast to existing global anomaly detectors, which require the entire cube. It does not use segmentation and it is less susceptible to spectral mixtures in local neighborhoods of the imagery.

The approach consists of three major parts: (i) a data dimension reduction method, which plays an important role on the overall approach, since it maps the data from their native multivariate space to a univariate domain in order to avoid the small sample size problem mentioned earlier, gaining in the process insensitiveness to illumination on objects, reducing the dominance of the blackbody response produced by Earth (if the longwave infrared region applies), among other benefits to be discussed later; (ii) a univariate semiparametric model [21], which is chosen because of its robustness to samples representing a mixture of material types; (iii) a parallel random sampling scheme that characterizes the unknown background clutter, using a binomial probability density function to model the likelihood of sampling targets by chance and erroneously labeling them as clutter, justifying multiple sampling trials in parallel in order to significantly decrease the labeling error probability.

The semiparametric model is neither parametric, since the specific distribution controlling the data is not assumed to be known, nor nonparametric, since other parameters must still be estimated relating two unknown distributions. The semiparametric model, however, assumes that the distributions of the samples to be tested are related to each other through an exponential function (a distortion), having two unknown parameters. As it will be shown later in detail, the model is appealing for many reasons, including the following: if two spectral samples under test belong to the same material type (i.e., they are not anomalous to each other), then the assumed exponential function relating both distributions is reduced to unity. If the two samples under test belong to different material types (i.e., they are anomalous to each other), then the exponential function will impose a significant weight relating both distributions, indicating that the samples are anomalous to each other; a key point here is that such an outcome is invariant to whether the assumption of exponential relationship between the distributions is satisfied or not—this will be discussed in more detail later. Finally, another benefit, although not recognized earlier by users of semiparametric models in other areas of study, is the model’s natural ability to handle samples representing a mixture of different material types, which also will be shown later. As a note, samples representing a mixture of material types are known to significantly increase the false alarm rate in operational scenarios, requiring autonomous anomaly detection; they can produce dominant edges between spatial regions of different material classes, later to be detected as meaningless (false) anomalies [2].

The strength of the semiparametric model handling the mixture problem is attributed to the fact that a sample under test is expected to contribute to the estimation of the distribution function controlling the sample labeled as reference, where both samples are expected to equally contribute to the estimation (when both are in fact under the same distribution), only the reference sample will be able to contribute (when both samples are from different distributions), or both the reference sample and a portion of the test sample will contribute (when the reference is a mixture and the test sample represents a component in that mixture, or vice versa). These outcomes are produced naturally by the model because of the imposed exponential relationship between the two distributions, as shown later using simulated univariate data to make the point. Another appeal for using an exponential distortion assumption is that many of the known distribution functions can be expressed in terms of an exponential distortion of another distribution, including all of the exponential family of distributions [22].

Rosario in [23] published a much earlier and limited conference paper version of this work, where a two-step semiparametric detector (data transformation and semiparametric test statistic) was introduced to the limited task of local anomaly detection (where prior knowledge of targets’ scales was required, imposing the limitation) and its performance was compared only against the RX algorithm. Relative to the previous work in [23], this paper significantly extends both the overarching methodology and presents additional results using the extended approach to test significantly more challenging scenarios from the ground to ground viewing perspective, where targets’ scales are unknown a priori. In other words, this work enables capability rather than just offers an incremental improvement. The specific contributions in this paper (method extension and additional results) are as follows. (i)The two-step anomaly detector (data transformation and semiparametric test statistic) is employed for the first time in a quasi-global framework (which is also proposed herein), where only a fraction of the entire data cube being represented by blocks of data are randomly selected from the imagery and used as reference sets of spectra to test the entire imagery. The results are later fused using order statistics, as the sampling scheme is modeled by the Binomial distribution.(ii)An analytical cutoff threshold is derived from the approach’s asymptotic cumulative probability of rejecting a null hypothesis, when either the null or the alternative hypothesis is true (given that the null hypothesis is based on a multisample testing and order statistic scheme).(iii)The extended approach is applied to additional real HS imagery to automatically find manmade objects in the scene, producing excellent results in difficult natural clutter scenarios viewed from both nadir and forward looking viewing perspectives.(iv)Performance of the extended approach is compared via ROC curves against multiple methods found in the literature, for example, global methods (k-means, Gaussian mixture model, global RX), local methods (kernel RX, standard RX, Fisher’s linear discriminant, and the local semiparametric detector discussed in [23]).(v)An analytical study of the two-sample test framework, using the local semiparameter detector in [23] as the base detector. (vi)An analytical study of the multi-sample fusion test framework, using the semiparametric model embedded in the quasi-global framework, which relaxes the prior knowledge requirement of target scales, hence, enabling scene anomaly detection tasks independently of the viewing perspective (nadir or forward looking).(vii)A study of the semiparametric model’s behavior in the presence of samples representing a mixture of two different material classes, which is the most common mixture case scenario given the sliding window sizes typically used in anomaly detection operational scenarios. (viii)A subsection specially devoted to discuss the motivation and give a qualitative assessment of the data transformation used in the two-step semiparametric anomaly detection of [23].

For convenience, a list of notations is available after the appendix.

This paper is organized as follows: Section 2 describes spatial data window models for HS sensing, a semiparametric model, and a hypothesis test; Section 3 discusses the sampling method, its probabilistic model, and introduces QG-SemiP; Section 4 discusses performance of QG-SemiP testing nadir and forward looking HS imagery, consisting of manmade objects in difficult natural background scenarios; Section 5 concludes the paper.

2. Problem Formulation, Data Transformation, and Semiparametric Model Description

The main goal of anomaly detection algorithms testing incoming imagery is to detect objects that are spectrally anomalous to its surroundings, yielding in the process a tolerable number of nuisance detections. In many surveillance and reconnaissance applications, it is desired that manmade objects are detected as being anomalous to the surrounding natural clutter. Both format and model of the data play a significant role in attempting to achieve the intended goal.

2.1. Remote Sensing Data and Data Format

Experiment was carried out on data sets from two distinct sensors and viewing perspectives: (i) the hyperspectral digital imagery collection experiment (HYDICE) sensor, from a nadir looking perspective; (ii) the SOC-700 hyperspectral sensor, from a forward looking perspective. Data from these sensors will be referred to in this paper as nadir looking imagery and forward looking imagery, respectively.

The HYDICE sensor records 210 spectral bands in the visible to near infrared (VNIR: bands between 0.38 and 0.97 μm) and short-wave infrared (SWIR: bands between 1.0 and 2.50 μm). An extended description of this dataset can be found, for example, in Schweizer and Moura (2000). The results shown in this section for one data subcube are representative for other sub-cubes in the HYDICE (forest radiance) dataset. An illustrative subcube (shown as an average of 150 bands; 640×100 pixels) is depicted in Figure 1 (Cube 1, top). We discarded water absorption and low signal to noise ratio bands; the bands used are the 23rd–101st, 109th–136th, and 152nd–194th. The scene consists of 14 stationary motor vehicles (targets near a treeline) in the presence of natural background clutter (e.g., trees, dirt roads, grasses). Each target consists of about 7×4 pixels, and each pixel corresponds to an area of about 1.3×1.3 square meters at the given altitude.

Figure 1 

Examples of nadir looking imagery (Cube 1) and forward looking imagery (Cubes 2 through 4). An effective anomaly detection algorithm suite would allow a machine to automatically detect the presence of manmade objects (targets), while suppressing the cluttered environment, using no prior information about what constitutes clutter background or target in the imagery.

The forward looking imagery used for this work was recorded using the SOC-700 VNIR HS spectral imager, which is commercially available off the shelf. The system produces HS data cubes of fixed dimensions 𝑅=640 by 𝐶=640 pixels by 𝐾=120 spectral bands between 0.38 and 0.97 μm. Figure 1 (Cube 2 through 4, bottom) depicts examples of the forward looking imagery, where each pixel in any of the three cube examples corresponds to the average of 120 band values. Data cubes Cube 2 and Cube 3 were collected during the summer of 2004 in California, USA; Cube 4 was collected during the spring of 2008 in New Jersey, USA. From actual ground truth, it is known that the scene in Cube 2 (see Figure 1) contains three motor vehicles and a standing person in the center of the scene (i.e., two pick-up trucks to the left in proximity to each other, a man slightly forward from the vehicles in the center, and a sport utility to the right). Although the natural clutter in Cube 2 and Cube 3 is dominated by Californian valley-type trees and/or terrain at the same general geolocation, the data in Cube 3 depict a significantly more complex scenario. From actual ground truth, it is known that in Cube 3 a sport utility vehicle is inconspicuously deployed in the shades of a large cluster of trees. Portions of the shadowed vehicle can be observed near the center in Cube 3. Cube 4 was recorded in a wooded region at Picatinny, where (according to the available ground truth) a sport car is located behind multiple tree trunks, partially obscuring the vehicle; see Figure 1 (Cube 4).

The four data cubes in Figure 1 are independently displayed as intensity images after linear mapping the gray scale of each to the range 0–255. Pixel intensities shown in each individual surface is only relative to corresponding values in that surface; in other words, pixel values representing the same material (e.g., general terrain) may be displayed with different intensities in another surface. This fact explains, for instance, the difference in brightness between the same terrain under similar atmospheric conditions shown in Cube 2 and Cube 3. (The strong reflections from certain parts of the vehicles captured by the sensor in Cube 2 are not as dominant in Cube 3 because the vehicle in Cube 3 is in tree shades; the open field in Cube 3 is then the strongest reflector in the scene).

Next, we present a model of observed data using a sliding 𝑛×𝑛 window in 𝐗 (a data cube). The data format of 𝐗 is shown in (1), where 𝑟(𝑟=1,…,𝑅) and 𝑐(𝑐=1,…,𝐶) index pixels 𝐱_𝑟𝑐 in the 𝑅×𝐶 spatial area 𝐗, where 𝑛≪𝑅 and 𝑛≪𝐶. Pixels within a fixed 𝑛×𝑛 block of data in 𝐗 are indexed from the upper left corner of this block using 𝑖𝑗 relative to rows and columns in 𝐗, where 𝑖=1,…,(𝑅−𝑛−1) and 𝑗=1,…,(𝐶−𝑛−1). A representation of an 𝑛×𝑛 window at pixel location (𝑖,𝑗)=(2,2) in 𝐗 is as follows:𝐱𝐗=¹¹,𝐱₁₂,…,𝐱_1𝐶𝐱₂₁,𝐱₍₂₊₁₎₁,⋮𝐱_{(2+𝑛−1)1},^{𝑛×𝑛window,whereinthiscase𝑖=𝑗=2}𝐱_𝑖𝑗,𝐱_𝑖(𝑗+1),…,𝐱_{𝑖(𝑗+𝑛−1)}𝐱_(𝑖+1)𝑗,𝐱_{(𝑖+1)(𝑗+1)},…,𝐱_{(𝑖+1)(𝑗+𝑛−1)}⋮𝐱_{(𝑖+𝑛−1)𝑗},𝐱_{(𝑖+𝑛−1)(𝑗+1)},…,𝐱_{(𝑖+𝑛−1)(𝑗+𝑛−1)},…,𝐱_2𝐶,…,𝐱_(2+1)𝐶⋮⋮,…,𝐱_{(2+𝑛−1)𝐶}𝐱⋮⋮⋮_𝑅1,𝐱_𝑅2,…,𝐱_𝑅C. (1)

Before pixels within a block of data can be used by a detector, they need to be rearranged to a sequence of multivariate variables. The rearrangement is made by concatenating individual rows in the 𝑛×𝑛 window in (1) as follows: 𝐖₁=𝐱_𝑖𝑗,…,𝐱_{𝑖(𝑗+𝑛−1)},𝐱_(𝑖+1)𝑗,…,𝐱_{(𝑖+1)(𝑗+𝑛−1)}𝐱,…,_{(𝑖+𝑛−1)𝑗},…,𝐱_{(𝑖+𝑛−1)(𝑗+𝑛−1)}=𝐲₁₁,…,𝐲_1𝑛1, (2) where 𝐖₁∈𝐑^𝐾×𝑛1, 𝑛₁=𝑛², and 𝐲_1ℎ∈𝐑^𝐾(ℎ=1,…,𝑛₁), such that 𝐲₁₁=𝐱_𝑖𝑗,𝐲₁₂=𝐱_𝑖(𝑗+1), and so forth until finally 𝐲_1𝑛1=𝐱_{(𝑖+𝑛−1)(𝑗+𝑛−1)}. Since a window can be anywhere in 𝐗 and 𝐗 represents any HS data cube, {𝐲_1ℎ}^𝑛1_ℎ=1 are considered random vectors and the entire set of spectra that constitutes 𝐗 will be observed through the 𝑛×𝑛 window.

Using the assumption that the random vectors in 𝐖₁ are independent and identically distributed (i.i.d.), the distribution of data within the window, using (2), can be simplified to the following: 𝐲₁₁,…,𝐲_1𝑛1∼i.i.d.𝑔₁(𝐲∣𝜽), (3) where 𝑔₁(𝐲∣𝜽) is a conditional multivariate probability density function (PDF) and 𝜽 is its parameter set; both 𝑔₁(𝐲∣𝜽) and 𝜽 are typically unknown for real applications.

Normally, an anomaly detector requires two input sets of spectra (𝐖₁∈𝐑^𝐾×𝑛1) and (𝐖₀∈𝐑^𝐾×𝑛0) to perform its task on 𝐗. The test sample (𝐖₁) is obtained at a fixed location 𝑖𝑗 in 𝐗, as shown in (1) and (2); but, the source to obtain a reference sample (𝐖₀) will depend on the application, or viewing perspective.

For the nadir looking imagery, the most popular sampling method is to use pixel vectors surrounding a 𝑛×𝑛 block of data to construct 𝐖₀, where 𝐖₁ is constructed using the block of data to be tested. Notice that this sampling method is not suitable for the forward looking imagery because a priori knowledge of target scales in the imagery are required to set the size of a separation (guard) region between the block of data to be tested and locally surrounding samples.

For the forward looking imagery, the reference input set 𝐖₀ could be made available from a spectral library, or be randomly selected from the testing data cube. In either case, 𝐖₀ would be a rearranged version of a 𝑛×𝑛 block of data. The latter is addressed in Section 3, where (in order to make such a test useful for real applications) 𝐖₁ is independently compared to multiple spectral sets 𝐖₀^(𝑓)∈𝐑^𝐾×𝑛0(𝑓=1,…,𝑁); fusing thereafter these comparison results in order to yield a decision (output) surface, as described in Section 3.

Both input sets 𝐖₁ and 𝐖₀ feed the anomaly detector.

As mentioned earlier, whether the viewing perspective is nadir or forward looking, mixtures of different materials in 𝐖₁ and/or in 𝐖₀ can significantly degrade anomaly detectors’ performances, as examined by Manolakis and Shaw in [2]. It is customary to assume normality in (3), or mixture of Gaussian distributions, but experience has shown (see [1]) that relaxation of these assumptions is needed.

We discuss next a two-step approach for anomaly detection, as introduced in [23], comprising of spectral transformation followed by the application of a univariate semiparametric model.

2.2. Data Transformation

This subsection discusses the employed method to transform spectra from their native multivariate space to a univariate domain, satisfying the univariate data requirement of the semiparametric model. We also provide justification for choosing the employed transformation and give some example cases to reinforce its use.

We consider a data transformation approach in two parts: (i) spectral differencing and (ii) angle mapping.

The rationale for (i) is twofold: (a) since HS samples are contiguous in the spectral domain (i.e., typical spectral resolution is of the order of 10 nanometers), more discriminant and independent information pertaining to a particular material type may be found between adjacent bands, which could augment the statistical power of detectors (this is specially the case in LWIR (longwave infrared) HS imagery where the radiance values observed in calibrated data (collected outdoor) are overwhelmingly influenced by the Planck’s blackbody equation [1]), as the Earth’s landscape (primarily, canopy and soil) behaves as a blackbody in the LWIR region of the electromagnetic spectrum (note: there is a whole topic of study in mathematical statistics on feature exploitation by zero crossing, which uses the output of random variable differencing as used herein for spectra, see [24]); (b) spectral differencing also puts significantly less weight in the importance of spectral magnitude (or bias) in anomaly detection applications, putting focus on the importance of spectral shape, instead. Spectral magnitude relates to the mean average of all measured radiance within a spectral sample, and spectral shape relates to the plotted curve of measured radiance as a function of wavelength. Existing classification and detection algorithms directly or indirectly exploit magnitude and/or shape of spectra in order to perform their tasks.

The benefit of (ii) is that it reduces the detection problem from a multivariate dimensional space to a univariate domain, avoiding the undesired problem of singularity issues during inverse estimations of covariance matrices. Singularity is known to occur when the sample size of a spectral sample is less than its number of spectral bands. Although there are approaches proposed in the literature to overcome this issue, the output of these approaches is not always desirable (see, for instance, [3]), since a typical HS sensor usually delivers between 120 and 1,000 bands, but targets may vary in number of pixels from as large as in the thousands to as small as 1 to 4 pixels, depending on the actual physical sizes of these targets and/or distance between the sensor and targets.

This paper is concerned about ensuring that the data transformation method can in fact reduce algorithm sensitivity to spectral magnitude (which can be achieved via (i)), so that a manmade object, for instance, deployed in tree shades can be considered as much as an anomaly relative to a dominant natural clutter background in the same way that the manmade object would have been if it were deployed, instead, out in the open field. All of this, while simultaneously preserving both a high sensitivity to spectral shape and the discriminant characteristics among spectra of distinct material types (which can be achieved via (ii)). If these requirements are matched, then the data transformation approach has achieved the overall desired goal—some example results are shown later in this subsection to reinforce those points.

Before those examples are presented, however, consider the following: let two spectra—having 𝐾 spectral bands—be available for comparison, 𝐲₀=[𝐿₀₁,…,𝐿_0𝐾] and 𝐲₁=[𝐿₁₁,…,𝐿_1𝐾], where 𝐿_𝑖𝑗(𝑖=0,1;𝑗=1,…,𝐾) are nonnegative radiance values. Further assume that 𝐲₁ is twenty percent stronger in magnitude than 𝐲₀. One way to formalize the disparity in magnitude is to let {𝐿_0𝑗=𝜇+𝛿_0𝑗}^𝐾_𝑗=1, {𝐿_1𝑗=1.2𝜇+𝛿_1𝑗}^𝐾_𝑗=1, and 𝜇>0. Two key points are worth noticing: (a) the difference 𝐿_0(𝑗+1)−𝐿_0𝑗=(𝜇+𝛿_0(𝑗+1))−(𝜇+𝛿_0𝑗)=𝛿_0(𝑗+1)−𝛿_0𝑗 would provide more discriminant and independent information (𝛿_0(𝑗+1)−𝛿_0𝑗) than jointly using the highly correlated 𝐿_0(𝑗+1) and 𝐿_0𝑗, ditto for 𝐿_1(𝑗+1)−𝐿_1𝑗; (b) the difference 𝐿_1(𝑗+1)−𝐿_1𝑗=𝛿_1(𝑗+1)−𝛿_1𝑗 would remove from consideration the 20 percent stronger average magnitude of 𝐿_1𝑗 over 𝐿_0𝑗, if the {𝐿_0(𝑗+1)−𝐿_0𝑗}^𝐾−1_𝑗=1 were used instead for comparison against {𝐿_1(𝑗+1)−𝐿_1𝑗}^𝐾−1_𝑗=1. In essence, (a) replaces the need for using—for instance—PCA (principal component analysis) to uncorrelate the data, as it is commonly employed in the HS community [2]); from (b), if 𝐲₁ and 𝐲₀ are observations of the same material under different illumination conditions (e.g., spectrum 𝐲₀ representing a shaded object and spectrum 𝐲₁ representing the same object but not shaded), then the average magnitude difference between the two spectra would not play a role in the comparison test, which is desired. For those readers who may have some concerns about what may be lost in the process of transforming the data from multivariate to univariate in the context of anomaly detection, as it will be shown shortly, the loss—although difficult to quantify—is not relevant to the anomaly detection problem, since an effective anomaly detector is not expected to distinguish a material type that is spectrally similar to another material type. If the detector is designed to be that sensitive, it would likely also produce an unacceptably high false alarm rate due to expected within class variability of the same material types dominating the scene (for instance, the expected within class variability of tree clusters across the scene).

The two-part transformation approach is described next.

Borrowing from the discussion in Section 2.1, the transformation approach requires two sets of spectra, a reference set (𝐖₀), 𝐖₀=𝐲₀₁,…,𝐲_0𝑛0=𝐿₀₁₁,…,𝐿_01𝑛0⋮𝐿_0𝐾1,…,𝐿_0𝐾𝑛0, (4) where 𝐲_0𝑖=[𝐿_01𝑖,…,𝐿_0𝐾𝑖]^𝑡 is calibrated spectra from a pixel-size location at the scene observed by the 𝐾-band sensor, during a particular set of atmospheric and illumination conditions; (⋅)^𝑡 is the vector transposed operator; 𝐿_0𝑘𝑖(𝑘=1,…,𝐾) are radiance values, such that, adjacent radiance values are usually highly correlated; 𝑖=1,…,𝑛₀ and 𝑛₀ is the sample size of 𝐖₀; and an independent test set (𝐖₁) that most likely has the same atmospheric condition captured in (4), but not necessarily the same illumination condition captured in (4), 𝐖₁=𝐲₁₁,…,𝐲_1𝑛1=𝐿₁₁₁,…,𝐿_11𝑛1⋮𝐿_1𝐾1,…,𝐿_1𝐾𝑛1, (5) where all of the qualifying comments made for (4) also apply to 𝐲_1𝑖=[𝐿_11𝑖,…,𝐿_1𝐾𝑖]^𝑡 with 𝑖=1,…,𝑛₁.

Letting 𝑢 denote the index that distinguish both sets 𝐖_𝑢(𝑢=0,1), the magnitude of 𝐿_𝑢𝑘𝑖 in (4) and (5) depends on the amount of illumination (e.g., shaded or nonshaded objects) and the illumination environment, this dependence can be significantly reduced by applying the first order differentiation—an approximation of the first derivative—to the columns of 𝐖_𝑢, or ∇₀=𝐿₀₂₁−𝐿₀₁₁,…,𝐿_02𝑛0−𝐿_01𝑛0⋮𝐿_0𝐾1−𝐿_0()1𝐾−1,…,𝐿_0𝐾𝑛0−𝐿_{0()𝑛𝐾−10},∇₁=𝐿₁₂₁−𝐿₁₁₁,…,𝐿_12𝑛0−𝐿_11𝑛1⋮𝐿_1𝐾1−𝐿_1(𝐾−1)1,…,𝐿_1𝐾𝑛0−𝐿_{1()𝑛𝐾−11}. (6)

Notice in (6) that ∇₀∈𝐑^{(𝐾−1)×𝑛0} and ∇₁∈𝐑^{(𝐾−1)×𝑛1}, and the sample means of ∇₀ and ∇₁ are, respectively, ∇₀=1𝑛₀∇₀𝟏_𝑛0×1,∇₁=1𝑛₁∇₁𝟏_𝑛1×1, (7) where 𝟏_𝑑×1 is a column vector of dimension 𝑑 filled with real values of 1’s.

Denoting the columns of ∇₀ (which corresponds to the reference set) as {∇_0𝑖∈𝐑^(𝐾−1)}^𝑛0_𝑖=1, then the multivariate reference and test samples can be transformed to univariate reference and test samples through the following angle-mapping formulas: 𝑥_0𝑖=180𝜋∇arccos^𝑡0𝑖∇₀∇_0𝑖Δ₀,𝑥_1𝑖=180𝜋∇arccos^𝑡0𝑖∇₁∇_0𝑖Δ₁, (8)(9) where 0^∘≤𝑥_0𝑖≤90^∘, 0^∘≤𝑥_1𝑖≤90^∘, the operator ‖⋅‖ using a column vector 𝐱 is the square root of 𝐱^𝑡𝐱 (note: although we prefer to use a metric that yields a number having a geometric interpretation, the reader who is concerned about algorithm speed may replace the angle mapper metric with any other comparable metric of choice, for instance, the correlation metric [2] or the normalized dot product showed inside the arccos operator in (8) and (9). The most important aspect about the employed metric is that it must be able to preserve the discriminant characteristics among spectra of different material types, as it will be shown shortly).

From (8) and (9), the following two univariate sequences are constructed: 𝑥₀=𝑥₀₁,𝑥₀₂,…,𝑥_0𝑛0,𝑥₁=𝑥₁₁,𝑥₁₂,…,𝑥_1𝑛0, (10) where 𝑥₀ (reference) and 𝑥₁ (test) are the input sequences to be used by the univariate based anomaly detection technique discussed in Section 2.3. Note that both samples end up having the same sample size, 𝑛₀.

As mentioned earlier, the employed data transformation was specifically chosen to offer a number of desired properties, including reduced sensitivity to spectral magnitude and preservation of discriminant features among spectra of distinct material types. Regarding the latter, notice in (10) that, for the proposed transformation to work as advertised, when both multivariate samples 𝐖₀ and 𝐖₁ happen to be observations of the same material type, the component values in 𝑥₀ and 𝑥₁ are expected to be comparable and closer to 0° in the scale between 0° and 90°; however, when 𝐖₀ and 𝐖₁ are observations of distinct material types, the component values in 𝑥₀ and 𝑥₁ should be proportionally apart, where values in 𝑥₀ are expected to be closer to 0° while values in 𝑥₁ are closer to 90°. In addition, when the observation represents a mixture of two different material types, the proposed transformation should yield a univariate sample that is representative of the mixture.

Now, we will take a closer look at these expectations, using (8) and (9) to transform two sets of real HS spectra for a qualitative assessment, addressing first the most common sliding window-based testing scenarios naturally occurring in anomaly detection problems: local testing, which requires a priori knowledge of object scales (range dependent), and global testing, which does not require a priori knowledge of object scales (range invariant). Figure 2 depicts these scenarios (in this context, local testing means comparing clustered spectra against neighboring spectra, while global testing means comparing clustered spectra against spectra located elsewhere across the same imagery).

Figure 2 

Most common window-based testing scenarios in anomaly detection problems, assuming for simplification that the scene background consists of two distinct material types (Class A and Class B) and a third material (Class C) also distinct from Class A and Class B depicts a anomalous material relative to the background.

Figure 2 illustrates the same data-cube spatial representation under the two-test case scenarios, where for simplification the scene background is spatially dominated by only two distinct material types (Class A and Class B) and a third material (Class C—also distinct from Class A and Class B) illustrates the presence of an anomalous material relative to the background. Notice also that the two objects of Class C in the scene have significant size and shape differences, so that, the advantage of approaching the anomaly detection problem from a global rather than a local perspective can be pointed out.

The left hand side image in Figure 2 shows the overlaid sliding window locations of the standard approach to local anomaly detection in the HS research community (see [3]), where a sliding window consisting of an interior square (inner window) is concentrically located along with a larger square so that spectra observed through the inner window can be compared against spectra observed through the outer portion (outer window) of the larger square (i.e., the area of the larger square minus the inner window). Furthermore, the spectral set observed through the outer window is labeled as the reference sample, while the spectral set observed through the inner window is labeled as the test sample.

As the test window slides across the image one pixel location per algorithmic test, the labels P1 through P7 (left hand side image in Figure 2) highlight seven key test locations in the image. Table 1 summarizes a list of plausible anomaly declarations by an anomaly detector versus the desired declarations, using the known ground truth information about the scene.

Table 1 

Plausible versus desired declarations of local anomalies.

The last column of Table 1 (desired anomaly) shows that out of the seven most distinct window locations for local testing, only two (P6 and P7) are desired, which may contradict declarations made by anomaly detection models currently found in the literature. In fairness, these detectors would be performing the job they are employed to do, as the plausible anomaly declarations shown in Table 1 are indeed correct in the strict sense. For instance, P2 shows a test between observations from Class A (test sample) and observations from a mixture (reference sample: Class A and Class B), while P6 shows a test between observations from the same material type (the intended anomaly: Class C). The reason, however, these declarations are not desired is that locations similar to P2 will likely accentuate edges as anomalies across the image, increasing the probability of false alarms, and locations similar to P6 will suppress the intended anomaly, decreasing the probability of correct detections. Location P6 also emphasizes the lack of robustness using the local testing approach to find anomalies in the scene, as a priori knowledge of object scales (consisting of the anomalous material) are not always available.

Regarding global testing shown in Figure 2 (the illustration at the right hand side), the window locations denoted as P1 through P5 depict distinct testing locations (observed test samples), while locations denoted as R1 through R3 depict distinct observations (fixed reference samples) that may have been randomly sampled from the image (prior to testing) and used to test every possible window-sized regions across the image, including P1 through P5. Table 2 shows the plausible declarations of anomalies versus the desired declarations, using the same ground truth information about the scene.

Table 2 

Plausible versus desired declarations of global anomalies (Using R1, R2, and R3 as reference samples, see Figure 2).

The last column of Table 2 (desired anomaly) in essence shows that any test that involves a mixture of classes and a component of that mixture should not be declared as an anomaly so that only truly anomalous material types (in this case Class C) relative to the background will be accentuated. An implementation scheme for the global testing approach, which requires fusion of declarations, will be discussed later. For now consider the following: the final declaration for any given window location is to retain the declaration NO, if it is there, out of all of the results produced by the particular testing location. Using this rule, locations P1, P2, and P3 would produce a final declaration of NO, and P4 and P5 would produce a final declaration of YES, as it would be desired by a global detection scheme. Notice also that P4 ensures that the circular object (Class C) would be accentuated as an anomaly, using the same test window size that would have detected the other anomaly of different scale shown in location P5; this is also desired.

Using real spectral samples, we will now qualitatively demonstrate the behavior of transformation output equations (10), when exposed to the key window positions shown in Figure 2. In anticipation, we would like to see that the data transformation will preserve distinction between samples of two different classes and produce results corresponding to a mixture of classes, when such a mixture is being observed. If successful, those results would give us some level of confidence that not much is being lost in terms of class distinction and that examples of mixtures would be shown as mixtures by the data transformation (demonstration of the desired anomaly detection declarations will be shown later, when the semiparametric model is discussed in Section 2.3. We also defer answering questions about what would happen when other possible window location cases appear, besides the ones shown in Figure 2, until discussion of results from testing real HS imagery in Section 4).

Figure 3 shows spectral transformation results using two sets of real spectra (Class A and Class B), where in this case spectral band differencing was not used as input for angle estimation among spectra and spectral means (see (6)), the actual individual radiance values (L’s) were used instead of the difference between adjacent radiance values (L’s). Figure 4 shows results using the same two sets of spectra but this time around using band differencing, following exactly the path from (6) to (10).

(a)

(b)

(c)

(d)

(e)

(f)

(g)

Figure 3 

Spectral transformation from multivariate to univariate sample using the actual spectral values instead of spectral differencing.

(a)

(b)

(c)

(d)

Figure 4 

Spectral transformation from a multivariate to univariate sample using spectral differencing.

Class A in Figure 3 consists of 200 real spectra representing a grassy area in an open field; Class B consists of 200 real spectra representing a cluster of tree leaves in the same geographic location of Class A. The employed HS sensor operates in the VNIR (visible to near infrared) region, producing 120 spectral bands per spectrum. To observe the behavior of the proposed data transformation as it processes spectra equivalent to the window location examples shown in Figure 2, Class A is denoted as the reference sample and processed with an independent test set representing Class A (200 spectra), it was also processed with a second test set representing Class B (200 spectra). The spectral sample means of both spectral sets are also shown in Figure 3.

Using both the reference and test samples from the same class (Class A) exemplifies window locations P1, P3, P5, and P6 in the local testing image in Figure 2; for the global testing illustration in Figure 2, it also represents the cases when spectra from R1 are processed with spectra from P1 and when spectra from R3 are processed with spectra from P3. Using spectra Class A as the reference and spectra from Class B as the test exemplifies the local testing locations P7 (i.e., the presence of an anomaly) and global testing location duals (R1, P4), (R1, P5), (R3, P4), and (R3, P5). The key point to notice in Figure 3 is the angle mapped plots shown in the bottom portion of the figure, where the plot at the left hand side shows that both the univariate reference sample (blue bubbles) and the univariate test sample (green bubbles) are comparable, preserving the fact that both samples belong to the same class. On the other hand, the angle mapped plot shown on the right hand side shows that the spectral transformation preserves the distinction between the samples from Class A and Class B. Both results are desired and would be passed as input to the employed detector for a decision.

The same experiment was held, but instead differentiating data in the spectral direction in order to check whether such a step would change the desired outcome from the spectral transformation shown in Figure 3. Figures 4(a) and 4(b) shows the band differencing means using (7) choosing samples from Class A to be both the reference and test sets (Figure 4(a)), and choosing a sample from Class A as reference and a sample of Class B as the test (Figure 4(b)).

The angle plots shown in Figures 4 (c) and 4(d) affirm that the differentiation step shown in (6) and (7), followed by output transformation from multivariate to univariate data, do preserve the lack of distinction between samples from the same class (Figure 4(c)) and a strong distinction between samples from Class A (reference) and samples from Class B (test), see plots in (Figure 4(d)). This example provides a direct assertion that what is lost due to the transformation is not relevant to the anomaly detection problem, since the goal of an anomaly detector is to find those material types that are truly distinct from the material types spatially dominating the background scene, as other factors will conspire against the detector’s effectiveness (e.g., mixture of distinct material classes, within material class variability). In other words, if a specific material type (target) is desired by the user to be considered as an anomaly relative to a natural environment but the target is not sufficiently distinct spectrally from the background, then an effective anomaly detector would most likely not declare the target as an anomaly. Ironically, as mentioned earlier, this is a virtue because a detector that is too sensitive in distinguishing two spectrally similar material types is also likely to produce an unacceptably high number of false alarms due to expected within class variability of dominant background material types in the scene.

Next, we would like to check how a spectral set, representing a mixture, is altered by the data transformation. For this demonstration, we constructed a reference sample by combining 100 spectra of Class A with 100 spectra of Class B, so that, the reference sample represented a mixture of two classes, and arbitrarily chose the test set to be represented by 200 spectra of Class B (the latter were independent from the ones used to construct the reference mixture). Figures 5(a) and 5(b) show both the mixture (reference construct) and the component of that mixture (Class B: test sample), and the resulting angle mapped plots using spectra without band differencing (Figure 5(c)) and spectra after band differencing (Figure 5(d)).

(a)

(b)

(c)

(d)

Figure 5 

Spectral transformation from a multivariate to univariate sample for a mixture and a component of that mixture.

The key message from Figures 5(c) and 5(d) is that the data transformation with or without the band differencing step does preserve in the univariate domain the fact that the material class of the multivariate test spectra is a component of the mixture of classes represented by the multivariate reference set of spectra. In this particular case, the univariate reference sample (blue bubbles) clearly shows the presence of two classes (i.e., half of the observations has lower angle values and the other half has higher angle values), while the test univariate sample (green bubbles) shows observation values commensurate to the one of the mixture class component (lower angle values). Although this fact is not necessarily desired for anomaly detection (see, for instance, P2 and P4 in Figure 2 (local testing, left image), the data transformation at least does not seem to alter such a scenario involving a mixture, which is fine as long as the employed anomaly detector is designed to handle similar challenging cases.

In summary, the proposed data transformation does preserve distinctions or similarities that exist between spectral sets and also offers additional benefits, as highlighted earlier in this subsection. What is needed now is a model that is more effective finding anomalies, while simultaneously managing the expected negative-impact nuisances naturally occurring in real operational scenarios (e.g., the presence of mixtures). A semiparametric model is described next for that purpose.

2.3. Univariate Semiparametric Model

Let two multivariate samples 𝐖₁ and 𝐖₀ be transformed to 𝑥₀=(𝑥₀₁,𝑥₀₂,…,𝑥_0𝑛0)∼𝑓₀(𝑥) and 𝑥₁=(𝑥₁₁,𝑥₁₂,…,𝑥_1𝑛1)∼𝑓₁(𝑥), respectively (using, for instance, (4) through (9)), where 𝑓₀(𝑥) and 𝑓₁(𝑥) are unknown joint PDFs.

To simplify the anomaly detection problem using the transformed data, suppose the two random samples (real valued, not vector valued) 𝑥₀ (reference) and 𝑥₁ (test) are independent, consisting of i.i.d. (see mathematical notation list after the appendix) random variables controlled by unknown marginal PDFs 𝑔₀ and 𝑔₁, respectively. Moreover, let 𝑔₀ and 𝑔₁ be related through the following model: 𝑥₁=𝑥₁₁,…,𝑥_1𝑛1i.i.d∼𝑔₁(𝑥)𝑥₀=𝑥₀₁,…,𝑥_0𝑛0i.i.d∼𝑔₀(𝑥),𝑔₁(𝑥)𝑔₀(𝑥)=exp𝛼+𝛽𝑥. (11)(12)

The model in (11)-(12) is appealing for many reasons, consider the following examples.

If 𝑥₀ and 𝑥₁ are samples of the same distribution, then the assumed exponential relationship is reduced to unity so that 𝑔₁=𝑔₀ (whether 𝑔₁ or 𝑔₀ is known or not), indicating that 𝑥₀ is not anomalous to 𝑥₁. If 𝑥₀ and 𝑥₁ are samples of different distributions, then the exponential function will impose a significant weight relating both distributions, indicating that 𝑥₁ and 𝑥₀ are anomalous to each other. The key point here is that the latter outcome is invariant to and independent of whether the assumption of exponential relationship between the distributions is satisfied or not; that is, if the exponential relationship assumption is satisfied, then 𝑥₁ is anomalous to 𝑥₀, but if this assumption is not satisfied, 𝑥₁ is still anomalous to 𝑥₀. Since the application of interest only requires a determination of whether 𝑥₁ is anomalous to 𝑥₀, satisfying the imposed assumption of (12) is inconsequential. So, a hypothesis test could just be designed to check whether exp(𝛼+𝛽𝑥)=1 in order to determine the presence of anomalies.

However, the relationship assumption in (12) plays a major role in the mathematical development of the statistical test and, more importantly for the application in context, it also plays an important role in determining whether a portion or the entire test sample (𝑥₁) can contribute to the estimation of the reference distribution (𝑔₀). The latter is a subtle feature never before recognized in other areas of study, for example, pharmaceutical [22], where semiparametric models are more commonly employed for their utility. The implication of this subtlety is that when one of the samples represents a mixture of different material types and the other represents a component of that mixture, the model in (11)-(12) allows both samples, through the assumed relationship, to contribute to the distribution estimation of the chosen reference sample; in essence stating that the assumption may also be partially satisfied as long as a portion of the test sample belongs to the reference distribution. This is manifested when the test produces an estimation of exp(𝛼+𝛽𝑥) that is significantly closer to unity than it would produce when the test sample has absolutely no relationship with the reference sample. In the practical scenario this is the difference between having a detector capable or not of naturally handling samples representing a mixture, as it will be shown later in this subsection. As mentioned earlier, samples representing a mixture of material types are known to significantly increase the false alarm rate in operational scenarios for autonomous anomaly detection—they can produce dominant edges between regions of different material classes, later to be detected as false alarms [2].

Notice in (12) that, since 𝑔₁ is a density, 𝛽=0 must imply that 𝛼=0, as 𝛼 merely functions as a normalizing parameter, following from the requirement that a PDF (in this case 𝑔₁) must integrate to unity, see PDF properties in [25]. Notice also that if 𝛽=0 then 𝑥₀ and 𝑥₁ must belong to the same population (i.e., 𝑔₁=𝑔₀). Using this fact, a test statistic can be designed to test the following hypotheses: 𝐻₀∶𝛽=0𝑔₁=𝑔₀𝐻anomalyabsent,₁∶𝛽≠0𝑔₁≠𝑔₀anomalypresent. (13)

In order to estimate 𝛽, denote the union of 𝑥₀ and 𝑥₁ (combined data) by 𝑡, 𝑡=𝑥₁₁,…,𝑥_1𝑛1,𝑥₀₁,…,𝑥_0𝑛0≡𝑡₁,…,𝑡_𝑛, (14) and following the construction by Qin and Zhang in [21] and Fokianos et al. in [22], a maximum likelihood estimator of 𝐺₀(𝑥)—the continuous cumulative distribution function (CDF) corresponding to the reference 𝑔₀(𝑥), can be obtained by maximizing the likelihood over the class of step CDF’s with jumps at the observed values 𝑡₁,…,𝑡_𝑛. Accordingly, if 𝑔₀(𝑡_𝑖)=𝑑𝐺(𝑡_𝑖), where 𝑑(⋅) denotes the differentiation operator, 𝑖=1,…,𝑛, the likelihood becomes, 𝜁𝛼,𝛽,𝑔₀=^𝑛0_𝑖=1𝑔₀𝑥_0𝑖^𝑛1_𝑗=1exp𝛼+𝛽𝑥_1𝑗𝑔₀𝑥_1𝑗=^{𝑛=𝑛1+𝑛0}_𝑖=1𝑔₀𝑡_𝑖^𝑛1_𝑗=1exp𝛼+𝛽𝑥_1𝑗. (15)

One can now express each 𝑔₀(𝑡_𝑖) in terms of 𝛼 and 𝛽 and then substitute the terms 𝑔₀(𝑡_𝑖) back into the likelihood to produce a function of 𝛼 and 𝛽 only. When 𝛼 and 𝛽 are fixed, (15) is maximized by maximizing only the product term ∏^𝑛_𝑖=1𝑔₀(𝑡_𝑖), subject to the constraints ^𝑛_𝑖=1𝑔₀𝑡_𝑖=1,^𝑛_𝑖=1𝑔₀𝑡_𝑖exp𝛼+𝛽𝑡_𝑖−1=0. (16)

Denoting 𝜌=𝑛₁/𝑛₀, Qin and Zhang in [21] showed that 𝑔₀𝑡_𝑖=1𝑛₀11+𝜌exp𝛼+𝛽𝑡_𝑖, (17) and that the value of the profile log-likelihood log[𝜁(𝛼,𝛽,g̃₀] up to a constant can be expressed as a function of 𝛼 and 𝛽 only, or 𝑙𝛼,𝛽=^𝑛1_𝑗=1𝛼+𝛽𝑥_1𝑖−^𝑛_𝑖=1log1+𝜌exp𝛼+𝛽𝑡_𝑖. (18)

A system of score equations that maximizes (18) over (𝛼,𝛽) is shown below [21], 𝜕𝑙𝛼,𝛽𝜕𝛼=−^𝑛_𝑖=1𝜌exp𝛼+𝛽𝑡_𝑖1+𝜌exp𝛼+𝛽𝑡_𝑖+𝑛₁=0,𝜕𝑙𝛼,𝛽𝜕𝛽=−^𝑛_𝑖=1𝑡_𝑖𝜌exp𝛼+𝛽𝑡_𝑖1+𝜌exp𝛼+𝛽𝑡_𝑖+^𝑛1_𝑗=1𝑥_1𝑗=0. (19)

The solution of the score equations yields the maximum likelihood estimators (𝛼,𝛽) and consequently by substitution also yields an estimator of 𝑔₀(𝑡_𝑖), or [21] 𝑔₀𝑡_𝑖=1𝑛₀11+𝜌exp𝛼+𝛽𝑡_𝑖. (20)

So, in summary, by using profiling, an estimator (20) for 𝑔₀(𝑡_𝑖) is attained in addition to score equations (19), where both the reference and test samples as shown in (14) are used to estimate 𝑔₀ (the reference PDF). This is only possible because the model in (12) implies that 𝑔₁ can be expressed in terms of 𝑔₀. This feature allows this model to be robust when either 𝑔₀ or 𝑔₁ is bimodal or multimodal (representing a sample mixture) while the other represents a component of the same mixture—a key factor for handling transitions of distinct regions in the anomaly detection application, as it will be shown later.

Using results from Fokianos et al. in [22], the estimator 𝛽 has the normal asymptotic behavior, as the sample size tends to infinity  (𝑛→∞), or √𝑛𝛽−𝛽₀→−_𝑛→∞𝑁𝜌0,⁻¹1+𝜌²𝑣², (21) where 𝛽₀ denotes the true parameter, 𝑣² is the variance (a scalor) using components from the combined sequence 𝑡, 𝜌=𝑛₁/𝑛₀, 𝑛=𝑛₁+𝑛₀, and →means converges to—in this case to a normal distribution having 0 mean and variance equals to 𝜌⁻¹(1+𝜌)²/𝑣².

Both estimators 𝛼 and 𝛽 are required to estimate 𝑣² in (21) via 𝑔₀(𝑡_𝑖). Denoting 𝑣² the estimator of 𝑣² and using results from [22], 𝑣²=_𝑖𝑡²_𝑖𝑔₀𝑡_𝑖−_𝑖𝑡_𝑖𝑔₀𝑡_𝑖², (22) where 𝑔₀(𝑡_𝑖) is shown in (20).

Normalizing the left side of (21) with 𝜌⁻¹(1+𝜌)²/𝑣², setting 𝛽_𝑜=0 and squaring the final result, and using known properties of the family of Chi square distributions [25], one can test 𝐻₀ in (13) with the following expression: 𝑍_SemiP=𝑛𝜌1+𝜌⁻²𝛽²𝑣²→−_𝑛→∞𝜒²₁, (23) which has the Chi square distribution asymptotic behavior with 1 degree of freedom, 𝜒²₁. Under the idealized assumptions of model (11) and (12), a decision can be based on the value of 𝑍_SemiP in (23), that is, high values of 𝑍_SemiP reject hypothesis 𝐻_𝑜, declaring the presence of local anomalies (note: Regarding (23), as typical from any asymptotic expression, the larger the value of 𝑛 is the more accurate and precise the approximation of the expression will be. Since 𝑛 in this context coincides with twice the sample size of the reference sample, sample sizes typically used in anomaly detection applications will suffice (greater than 100, yielding 𝑛 greater than 200). Practitioners in statistics usually recommend that for univariate variables, asymptotic expressions should use at least thirty two observations, indicating that 𝑛 in this case should be at least 32 or greater).

The test statistic in (23) will be referred to from here forward as the SemiP test statistic or SemiP anomaly detector, which has two steps: data transformation and test statistic estimation.

2.4. Implementation Notes for the Standing Alone 𝑆𝑒𝑚𝑖𝑃 Detector

2.4.1. Function Maximization

In order to implement (23), we perform an unconstrained maximization of the log maximum likelihood function in (18), or alternatively one could minimize the negative version of (18), to obtain the extremum estimators 𝛼 and 𝛽. We used one of the conventional unconstrained nonlinear optimization algorithms—the simplex method [26], which is available in Matlab software (Release 13) under the function name fminsearch. The simplex method is a direct search method that does not use numerical or analytic gradients. If 𝑛 is the length of 𝑥, a simplex in 𝑛-dimensional space is characterized by the 𝑛+1 distinct vectors that are its vertices. For instance, in two-space, a simplex is a triangle; in three-space, it is a pyramid. At each step of the search, a new point in or near the current simplex is generated. The function value at the new point is compared with the function’s values at the vertices of the simplex and, usually, one of the vertices is replaced by the new point, giving a new simplex. This step is repeated until the diameter of the simplex is less than the specified tolerance. A limitation using such a method is that it may find a local extremum, so the choices of initial parameters may prove to be critical in some cases; however, we found in practice that by setting the initial values to (𝛼=0,𝛽=0), the method converges reasonably fast and works very well for all of the cases that we have observed, independently of whether anomalies were present or not in the tests.

2.4.2. Variance under the Null Hypothesis

The term 𝑣² in (23) is computed using 𝑣²=𝐸(𝑡²)−𝐸²(𝑡), where 𝐸(𝑡^𝑘)=∑_𝑖𝑡^𝑘_𝑖𝑔₀(𝑡_𝑖) and 𝑔₀(𝑡_𝑖) is shown in (20).

2.4.3. Decision Threshold

As mentioned earlier, using (23), high values of Z_SemiP reject hypothesis 𝐻_𝑜 in (13), declaring that 𝑥₁ is an anomaly relative to 𝑥₀. Using this detector as a standing alone unit, one could set a decision threshold based on the type I error, that is, based on the probability of rejecting 𝐻_𝑜 given that 𝐻_𝑜 is true. Using a standard integral table for the Chi square distribution, with 1 degree of freedom, find a threshold that yields an acceptable probability of error (e.g., 0.001).

2.5. Model Behavior in the Presence of Sample Mixtures

We show in this subsection the robustness of the semiparametric model toward an asymmetric test, that is, when a sample of a mixture is compared against a component of that mixture, which is found locally across the image in the form of spatial transitions. More specifically, we would like to show that 𝛽 (estimator for 𝛽 in (13)) is significantly closer to ZERO when, for two PDFs 𝑔_𝐴(𝑦)≠𝑔_𝐵(𝑦), 𝑦₁∼𝑔_𝐵(𝑦) and 𝑦₀∼𝑔_𝐵(𝑦) or 𝑦₀∼(𝑔_𝐴(𝑦)∪𝑔_𝐵(𝑦)) (representing the union ∪ or a mixture of two PDFs) than when 𝑦₁∼𝑔_𝐵(𝑦) and 𝑦₀∼𝑔_𝐴(𝑦). We illustrate this fact using simulated data and focusing on three specific case studies.

Case 1. 𝑦₀∼𝑔_𝐴(𝑦) versus 𝑦₁∼𝑔_𝐵(𝑦).

Case 2. 𝑦₀∼(𝑔_𝐴(𝑦)∪𝑔_𝐵(𝑦)) versus 𝑦₁∼𝑔_𝐵(𝑦).

Case 3. 𝑦₀∼𝑔_𝐵(𝑦) versus 𝑦₁∼𝑔_𝐵(𝑦).

According to [20], Case 2, which represents a transition of distinct regions in the image, appears some 20% of the entire image, or higher, as local patches are observed through a small moving window across typical images. Therefore, Case 2 is a major source of false alarms that could be avoided using a more robust model of the background than the typical models used in the target community. In anomaly detection applications, it is desired that the employed detector declares the presence of an anomaly for Case 1 but no anomalies for Cases 2 and 3; a declaration of no anomalies present is also desired, if Case 2 were reversed, that is, 𝑦₀∼𝑔_𝐵(𝑦) versus 𝑦₁∼(𝑔_𝐴(𝑦)∪𝑔_𝐵(𝑦)), although this case is not shown here, its results are consistent with Table 3.

Table 3 

Behavior of 𝛽 on different case studies.

The results shown in Table 3 were computed using 100 simulated random samples from an i.i.d. Gaussian distribution to represent the reference sequence 𝑦₀ and another 100 samples to represent the test sequence 𝑦₁. A sequence representing a mixture of two classes consists of 50 samples for each class resulting in a total of 100 samples. The formulation and parameters used to generate these sequences are shown in Table 3 for different case studies, where the samples in row 2 (starting from the left upper corner in Table 3) simulates a local test between a genuine isolated object (𝑦₁) and its homogeneous surrounding background (𝑦₀)—Case  1, the samples in row 3 simulates a local test at a transition between two classes, where the test sample belongs to one of these classes—Case  2; the samples in row 4 simulates a local test within a homogeneous region—Case  3. Practical implementation details of the SemiP detector, which includes the estimation of (𝛼,𝛽), are shown in Section 2.3. The parameters (𝛼,𝛽) were estimated by maximizing the log likelihood function, using an optimization subroutine initialized to (0,0), (0,0), and (0,0) for Cases  1,  2, and 3, respectively, so that convergence to a solution down to a tolerable error could be achieved using the subroutine.

Since the solution of the semiparametric model uses the union of samples 𝑡 and estimators 𝛼 and 𝛽 to estimate 𝑔₀ (which itself is an estimator of 𝑔₀), we also included in Table 3 the mean estimates 𝜇, 𝜇₂, and 𝜇₃ and the variance estimates 𝑣², 𝑣²₂, 𝑣²₃; where 𝑣² estimates variance from the solution of the semiparametric model using the union of samples 𝑡=(𝑦₀,𝑦₁) and 𝑔₀; 𝑣²₂ estimates variance using only the reference sample 𝑦₀; and 𝑣²₃ is the sample variance using 𝑡. The mean estimates were computed accordingly, see Table 3.

In reference to results shown in Table 3, recall that the null hypothesis is 𝛽=0, and notice that the value of 𝛽 in Table 3 are significantly closer to zero for Case  3 (homogeneous region) and Case  2 (a transition of two different region) than for Case  1 (genuine local anomaly), where in Case  1 𝑦₀ and 𝑦₁ do belong to different classes. Notice also that the disparity between the values of 𝑣²₂ and 𝑣²₃ for each case study also reflects how close 𝛽 is to zero. For instance, the disparity between 𝑣²₂, and 𝑣²₃ for Case  1 is quite large compared to corresponding disparities for Cases  2 and 3.

The semiparametric model handles mixture by showing sensitivity on the estimation of 𝛽 and 𝛼, such that when the test sample has strong statistical information about one of the subclasses in the reference sample, the semiparametric method responds by keeping both 𝛽 and 𝛼 relatively close to zero in order to maximize the log likelihood function in (18). The estimates 𝛼,𝛽 affect directly the computation of 𝑔₀, which in turn is used to compute 𝑣².

To shed some light on the effect of 𝛽 and 𝛼 on 𝑔₀, we present some results in Figure 2. The plots shown in Figure 2 (row 1) corresponds to Case 1; where, the plots on the left depict the values of 𝑡_𝑖 as a function of index 𝑖, for convenience we have marked where the sequences 𝑦₀ and 𝑦₁ are relative to each other within 𝑡; and the plots on the right depict 𝑔₀(𝑡_𝑖) versus 𝑖. Likewise for Case 2 (row 2) and Case 3 (row 3).

Let us consider first Cases 3 and 1. As mentioned earlier, because of the semiparametric model in (11) and (12), the union of samples 𝑡=(𝑦₀,𝑦₁), where 𝑦₀ is the reference sequence, is used to estimate the reference PDF estimator 𝑔₀. Circumstances when both samples belong to the same population (Case 3), the estimated cumulative weight for the test sample 𝑦₁, that is, 𝑤₁=∑^𝑛1_𝑖=1𝑔₀(𝑦_1𝑖), is expected to be approximately equal to the estimated cumulative weight for the reference sample 𝑦₀, that is, 𝑤₀=∑^𝑛1+𝑛0_𝑖=𝑛1+1𝑔₀(𝑦_{0(𝑖−𝑛1)}); because the constraint ∑^{𝑛=𝑛1+𝑛0}_𝑖=1𝑔₀(𝑡_𝑖)=1 was used in the profiling method to attain 𝑔₀(𝑡_𝑖) in terms of 𝛼,𝛽, we would expect both 𝑤₁ and 𝑤₀ to be near 0.500. Using simulated samples for the normal distributions and parameters shown in Table 3 for Case 3, we obtained 𝑤₁=0.4998 and 𝑤₀=0.5002, which are very close to our expectations, see Figure 6. We interpret the actual values of 𝑔₀(𝑡_𝑖) shown for Case 3 in Figure 6(f) to be an indication that the semiparametric method regards the test sample 𝑦₁ to be as important in the estimation of 𝑔₀ as the reference sample 𝑦₀ is in that estimation, for the right justification, as both 𝑦₀ and 𝑦₁ happens to be governed by the same distribution.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6 

Effect of estimators 𝛼,𝛽 on the computation of 𝑔0.

In contrast, when both 𝑦₀ and 𝑦₁ belong to clearly different classes (see Case 1 in Table 3), the semiparametric method recognizes this difference and virtually shuts down the contribution of 𝑦₁ in the estimation of 𝑔₀. The way it shuts down the contribution of 𝑦₁ is by maximizing the log likelihood function with values of 𝛽 and 𝛼 that allow the estimates of atomic exponential distortions exp(𝛼+𝛽𝑡_𝑖) to be relatively high when 𝑡_𝑖 corresponds to components of 𝑦₁.

And since exp(𝛼+𝛽𝑡_𝑖) are inversely proportional to 𝑔₀(𝑡_𝑖), see (20), the contributions of corresponding components of 𝑦₁ in 𝑡_𝑖 estimating 𝑔₀ are shut down as an indication of nonimportance to this estimation. The implication of this shut down is that the value of 𝛽 is relatively away from zero (relative to Case 3, for instance), which rejects the null hypothesis as desired. Figures 6(a) and 6(b) show the combined samples for Case 1 and the resulting cumulative weights for the test sample 𝑤₁=∑^𝑛1_𝑖=1𝑔₀(𝑦_1𝑖)≅0.0 and for the reference sample 𝑤₀=∑^𝑛1+𝑛0_𝑖=𝑛1+1𝑔₀(𝑦_{0(𝑖−𝑛1)})≅1.0, where ≅ denotes approximately equal to. The shutdown is reflected in the results for 𝑤₁.

In reference to Case  2, where the information carried in 𝑦₁ is also contained in 𝑦₀, as half of the random components in 𝑦₀ are governed by the same distribution of the random components in 𝑦₁—see Table 3 and Figures 6(c) and 6(d), the semiparametric method recognizes this fact by holding the value of 𝛽 at near zero, and interestingly by allowing the contributions of 𝑦₁ estimating g̃₀ to be comparable to those contributions of the portion of 𝑦₀ that are similar to 𝑦₁. In other words, since both 𝑦₁ and 𝑦₀ are used to estimate 𝑔₀ and half of the random components in 𝑦₀ are governed by the same distribution of all of the random components in 𝑦₁ and the other half are governed by a different distribution, the semiparametric method will not discriminate between 𝑦₁ and the portion of 𝑦₀ that is similar to 𝑦₁. The outcome of this behavior is that, in order to maximize the log likelihood function, the values of 𝛽 and 𝛼 are kept statistically close to zero (in this case 𝛽=−0.0073) to reflect that the information in 𝑦₁ is important in the estimation of 𝑔₀. The way this method explicitly shows this nondiscrimination is by yielding the cumulative weight using 𝑦₁ and the portion of 𝑦₀ that is similar to 𝑦₁ to hold half of the power, while the other half of the power is allocated to the portion of 𝑦₀ that is dissimilar to 𝑦₁. In other words, using results from Figure 6(d) and this claim, we would expect a value of 0.5 (half power) for the cumulative weight using 𝑦₁ and the portion of 𝑦₀ that is similar to 𝑦₁; we obtained ∑^𝑛1_𝑖=1𝑔₀(𝑦_1𝑖)+∑^𝑛1+𝑛0_{𝑖=𝑛1+51}𝑔₀(𝑦_{0(𝑖−𝑛1)})=0.5007. And we would expect the other half of the power to be in the cumulative weight using the portion of 𝑦₀ that is dissimilar to 𝑦₁; we obtained ∑^{𝑛1+𝑛0+50}_𝑖=𝑛1+1𝑔₀(𝑦_{0(𝑖−𝑛1)})=0.4993. Notice that adding 0.5007 and 0.4993 yields exactly 1.0 as expected because the constraint ∑^{𝑛=𝑛1+𝑛0}_𝑖=1𝑔₀(𝑡_𝑖)=1 was used in the profiling method in order to attain a representation of 𝑔₀ in terms of free parameters 𝛽 and 𝛼.

A conclusion that we can draw from this discussion is that the semiparametric method will indirectly compare two samples 𝑦₀ (reference) and 𝑦₁ (test) by assuming that the distribution of 𝑦₁(𝑔₁) and the distribution of 𝑦₀(𝑔₀) are related (exponentially) to each other and that, therefore, the information content in both samples can be used to estimate one of these distributions, in particular, 𝑔₀.

We found this indirect comparison method to be highly sensitive to the cumulative contribution of 𝑦₁ estimating 𝑔₀. This sensitivity has an important practical value in the anomaly detection application for three reasons.

First, if 𝑔₁=𝑔₀, both samples 𝑦₀ and 𝑦₁ are expected to equally contribute to the estimation of 𝑔₀, which in fact would improve that estimation due to the increase of sample size. Result: 𝑦₁ would be labeled as not being anomalous to 𝑦₀ in this application.

Second, if 𝑔₁≠𝑔₀, sample 𝑦₁ is not expected to be allowed to contribute to the estimation of 𝑔₀, thus, this estimation would solely rely on the cumulative contribution of 𝑦₀. Result: 𝑦₁ would be labeled as being anomalous to 𝑦₀ in this application.

And third, if 𝑔₀ is a mixture of densities, such that, 𝑔₁ is a component in that mixture, we found that the contribution of 𝑦₁ would not be suppressed, but proportional to the weight of 𝑔₁ in that mixture (see Figure 6). Result: 𝑦₁ would be labeled as not being anomalous to 𝑦₀ in this application.

This behavior of the semiparametric test statistic is highly desired in the target community because it conforms with be behavior of a human analyst performing the same task in the target application, and it separates this method from existing ones performing the same task.

3. Quasiglobal Semiparametric Approach

The semiparametric test statistic is used as the primary scoring method for the quasi-global anomaly detection approach. As mentioned in Section 1, the quasi-global anomaly detection approach was designed to tackle the forward looking anomaly detection problem, although the application of the quasi-global algorithm using nadir looking imagery are also considered in Section 4.

We start by describing the background sampling method and its probabilistic model, followed by description of the quasi-global algorithm framework using the semiparametric test statistic.

3.1. Sampling Method and Its Probabilistic Model

Assume that target pixels are present in the 𝑅×𝐶 spatial area of a 𝑅×𝐶×𝐾 HS data cube 𝐗, denote 𝑎 the total number of target pixels in 𝐗, 𝑞 the probability of a pixel in 𝐗 belonging to the target, and the relationship 𝑞=𝑎/𝐴, where 𝐴=𝑅𝐶 (or all pixels in 𝐗) (in most applications 𝑞 is unknown, and if multiple targets are present in the imagery, 𝑎 will be the total number of all pixels belonging to all targets present in the imagery; also, these targets may or may not have the same material type). In order to represent the unknown clutter background in the imagery, let 𝑁 blocks of data—all having a fixed small area (𝑛×𝑛)≪(𝑅×𝐶)—be randomly selected from the 𝑅×𝐶 area, see one of the data cubes in Figure 1. In theory, for (𝑛×𝑛)=(1×1) and using the assumption that target pixels in 𝐗 are disjoint and randomly located across the 𝑅×𝐶 imagery area (note that in practice, this assumption is not satisfied when targets are present in the scene, but we will use this assumption to establish a guideline), the probability 𝑃 that at least one block of data has a target pixel is 𝑃(𝑚≥1)=1−𝑝(𝑚=0), (24) where 𝑝 is the binomial density function [27], given parameters 𝑞 and 𝑁, and 𝑚∈{0,1,…,𝑁} is the number of blocks of data containing a target pixel, or 𝑝𝑚∣𝑞,𝑁=𝑁!()!𝑞𝑚!𝑁−𝑚^𝑚1−𝑞^𝑁−𝑚, (25) (symbols | and ! denote given parameters and the factorial operator, resp.).

For convenience, we will refer to 𝑃(𝑚≥1) as the probability of contamination and 𝑚 the number of contaminated blocks of data.

The implementation of this contamination model to the autonomous background sampling problem requires that each one of the 𝑁(𝑛×𝑛) blocks of data be regarded as an independent reference set 𝐖₀^(𝑓)(𝑓=1,2,…,𝑁) representing clutter spectra, where 𝐖₀^(𝑓)∈𝐑^𝐾×𝑛0 is a rearranged sequence version of the 𝑓th block of data having 𝑛₀=𝑛² spectra. By necessity, 𝑛₀ must be significantly greater than one—for statistical purposes—but yet significantly smaller than 𝐴=𝑅𝐶 (e.g., 𝑛₀/𝐴=20²/640²=0.000977) in order to reasonably satisfy the assumption that a 𝑛×𝑛 block of data has an unit area at the 𝑅×𝐶 imagery area. A contaminated block of data, then, will be treated qualitatively as a block having target pixels covering a large portion of the block’s area (e.g., greater than 0.70). In addition—when targets are present, since pixels representing a single target are expected to be clustered in the imagery, the assumption that each target pixel is randomly located across the imagery area will be ignored. Using (24), while ignoring the nonclustered target pixel assumption, implies that the probability of contamination will be overestimated, as blocks of data are less likely to be randomly selected from the same cluster of target pixels (for the autonomous background sampling problem, it is more conservative to overestimate the probability of contamination than to underestimate).

Figure 7(b) shows a plot of the probability of contamination 𝑃(𝑚≥1) versus 𝑁, for two values of 𝑞 (0.1 and 0.2). It is highlighted in the plot in reference that, for instance, if parameters are set to (𝑞,𝑁)=(0.10,22) then𝑃(𝑚≥1)=0.90. Notice that for 𝑁=22, if target pixels are present but cover less than 𝑞=0.10 of the imagery area, 𝑃(𝑚≥1) = 0.90 is overestimated by two fronts: (i) pixels from a single target are not randomly spread across the imagery area, but clustered, and (ii) the cumulative number of target pixels covers less than 0.10 of the imagery area. So, (24) provides an upper bound (conservative) approximation of the probability of contamination, given parameters 𝑞 and 𝑁.

(a)

(b)

(c)

Figure 7 

The relationship between 𝑁(the number of randomly selected blocks of data, shown as yellow squares on the imagery) and the contamination probability 𝑃𝑔(𝑚≥1)=𝑃(𝑚≥1) is shown in (b) for a given 𝑞 (e.g., 𝑞=0.10), which is an upper bound guess representing the maximum ratio between target pixels over the 𝑅×𝐶 area. To better characterize the unknown clutter background, a high 𝑁is most desired, but at a high cost, that is, an undesirably high 𝑃𝑔(𝑚≥1). The overall contamination probability, however, can significantly decrease by independently repeating the random sampling process 𝑀number of times, as shown (c) of figure, and then fusing results using a suitable method.

Figure 7(b) also shows the tradeoff between having a larger number of spectral sets (increasing 𝑁) in order to adequately represent the clutter background, which is desired, and the cost of increasing probability of contamination, which is not desired (in particular, contamination implies that once target pixels are randomly selected by chance from the imagery area, they will be used by a detector as reference set to test the entire imagery, in which case targets would be suppressed).

Since the presence of target pixels in the imagery is unknown a priori, finding a way to decrease the probability of contamination becomes a necessity. In order to decrease this probability, using an adequately large 𝑁 and a sensible value for 𝑞, we propose to independently repeat the random sampling process described in this subsection 𝑀 number of times. Figure 7(c) illustrates the outcome of 𝑀 repetitions. If we denote the probability of contamination of the 𝑔th random sampling process (or repetition) as 𝑃_𝑔(𝑚≥1), 1≤𝑔≤𝑀, for a fixed 𝑞 and 𝑁, note that each 𝑃_𝑔(𝑚≥1)=𝑃(𝑚≥1) and, since 0.0≤𝑃(𝑚≥1)≤1.0 and these processes will be repeated independently from each other, the overall probability 𝑃 that all of the processes will be contaminated with at least a contaminated block of data will decrease as a function of increasing 𝑀, or 𝑃=𝑃₁(𝑚≥1)𝑃₂(𝑚≥1)⋯𝑃_𝑀(𝑚≥1)=[𝑃(𝑚≥1)]^𝑀. (26)

Figure 7(c) plots 𝑃 as a function of increasing 𝑀, for 𝑃(𝑚≥1)=0.90 and 𝑃(𝑚≥1)=0.65. Taking, as an example, the 𝑃 curve in Figure 7 corresponding to using 𝑃(𝑚≥1)=0.90 in (24), notice that for 𝑀>40, [𝑃(𝑚≥1)]^𝑀 decreases to virtually zero. This outcome implies that at least one out of the 𝑀>40 processes has an extremely high probability of not being contaminated, given that 𝑁=22 and target pixels do not cover significantly more than 10% of the imagery area (𝑞=0.10).

3.2. Algorithmic Fusion

A framework for the quasi-global semiparametric anomaly detection algorithm can now be developed using (i) the repeated random sampling model discussed in Section 3.1 (needed to characterize the unknown clutter background in the imagery), (ii) the semiparametric anomaly detector discussed in Section 2.3 (needed to test reference data against the entire imagery), (iii) a way to fuse the results from testing 𝑁 randomly chosen blocks of data against the entire imagery using small windows (this will produce a 2-dim output surface per repetition), and (iv) a way to fuse 𝑀 independently produced 2-dim output surfaces into a single 2-dim decision surface.

Start by letting a HS data (𝑅×𝐶×𝐾) cube 𝐗 be available for autonomous testing. Let also 𝑁 blocks (𝑛×𝑛) of data be randomly selected from the 𝐗’s 𝑅×𝐶 spatial area and used as a reference library set 𝐖₀^(𝑓)(𝑓=1,2,…,𝑁) representing clutter background spectra, where 𝐖₀^(𝑓)=(𝐲^(𝑓)₀₁,…,𝐲^(𝑓)_0𝑛0) is a rearranged sequence version of the 𝑓th block of data having 𝑛₀=𝑛² spectra, where{𝐲^(𝑓)_0𝑢}^𝑛0_𝑢=1∈𝐑^𝐾 are 𝐾-dim column vectors. Let 𝐖₁=(𝐲₁₁,…,𝐲_1𝑛1) be the rearranged version of a (𝑛×𝑛) window of test data at location 𝑖𝑗 in 𝐗—see (1) for column vectors {𝐲_1ℎ}^𝑛1_ℎ=1∈𝐑^𝐾; first, we would like to automatically test 𝐖₁ against all {𝐖₀^(𝑓)}^𝑁_𝑓=1, and produce a single output (scalar) value 𝑍^(𝑖𝑗)_SemiP≥0.0 from these 𝑁 test results.

For better clarity in this subsection, we repeat the data transformation steps discussed in Section 2.2, but with the inclusion of index 𝑓=1,…,𝑁 and letting 𝐲_𝑖𝑗=𝐿_𝑖1𝑗,…,𝐿_𝑖𝐾𝑗, where 𝐿_𝑖𝑘𝑗 is the 𝑘th radiance value in 𝐲_𝑖𝑗, 𝑘=1,…,𝐾, and 𝐖₀^(𝑓)=𝐲^(𝑓)₀₁,…,𝐲^(𝑓)_0𝑛0=𝐿^(𝑓)₀₁₁,…,𝐿^(𝑓)_01𝑛0⋮𝐿^(𝑓)_0𝐾1,…,𝐿^(𝑓)_0𝐾𝑛0,∇₀^(𝑓)=𝐿^(𝑓)₀₂₁−𝐿^(𝑓)₀₁₁,…,𝐿^(𝑓)_02𝑛0−𝐿^(𝑓)_01𝑛0⋮𝐿^(𝑓)_0𝐾1−𝐿₀₍₎₁^(𝑓)_𝐾−1,…,𝐿^(𝑓)_0𝐾𝑛0−𝐿_0()𝑛^(𝑓)_𝐾−10,∇₀^(𝑓)=1𝑛₀∇₀^(𝑓)𝟏_𝑛0×1 (27)(28)(29) and, denoting the columns of ∇₀^(𝑓) as {∇^(𝑓)_0𝑢}^𝑛0_𝑢=1, 𝑥^(𝑓)_0𝑢=180𝜋∇arccos^{(𝑓)0𝑢𝑡}∇₀^(𝑓)∇^(𝑓)_0𝑢Δ₀^(𝑓)𝑛0_𝑢=1. (30) And equivalently for 𝐖₁=(𝐲₁₁,…,𝐲_1𝑛1)—the rearranged version of a (𝑛×𝑛) window of test data at location 𝑖𝑗 in 𝐗 and the columns of ∇₀^(𝑓) in (28)—{∇^(𝑓)_0𝑢}^𝑛0_𝑢=1, one has 𝑥^(𝑓)_1𝑢=180𝜋∇arccos^{(𝑓)0𝑢𝑡}∇₁∇^(𝑓)_0𝑢Δ₁^𝑛0_𝑢=1. (31) From (30) and (31), the following two univariate sequences will be used as inputs to the SemiP detector: 𝑥₀^(𝑓)=𝑥^(𝑓)₀₁,𝑥^(𝑓)₀₂,…,𝑥^(𝑓)_0𝑛0𝑥_1(𝑓),=𝑥^(𝑓)₁₁,𝑥^(𝑓)₁₂,…,𝑥^(𝑓)_1𝑛0, (32)(33) where 1≤𝑓≤𝑁.

Following the discussion that led to (33), results from the semiparametric test statistic can be used (or fused) as following: 𝑍^(𝑖𝑗)_SemiP=min_{1≤𝑓≤𝑁}𝑍^{(𝑖𝑗)(𝑓)}_SemiP, (34) where 𝑍^{(𝑖𝑗)(𝑓)}_SemiP=𝑛𝜌1+𝜌⁻²𝛽^(𝑓)2𝑣^2(𝑓), (35){𝑍^{(𝑖𝑗)(𝑓)}_SemiP}^𝑁_𝑓=1≥0.0, 𝑛₁=𝑛₀=𝑛², 𝜌=𝑛₁/𝑛₀, (𝑖=1,…,𝑅−𝑛−1) and (𝑗=1,…,𝐶−𝑛−1) index the left-upper corner pixel of an 𝑛×𝑛 window in 𝐗; using (32) and (33), 𝑡^(𝑓)=𝑥^(𝑓)₀₁,…,𝑥^(𝑓)_0𝑛0,𝑥^(𝑓)₁₁,…,𝑥^(𝑓)_1𝑛0=𝑡₁^(𝑓),…,𝑡_𝑛^(𝑓)_1+𝑛0,𝑣^2(𝑓)=^{𝑛1+𝑛𝑜}_𝑢=1𝑡²_𝑢^(𝑓)𝑔₀𝑡_𝑢^(𝑓)−^𝑛1+𝑛0_𝑢=1𝑡_𝑢^(𝑓)𝑔₀𝑡_𝑢^(𝑓)2,𝑔₀𝑡_𝑢^(𝑓)=1𝑛₀11+𝜌exp𝛼^(𝑓)+𝛽^(𝑓)𝑡_𝑢^(𝑓), (36)(37)(38) and estimates 𝛼^(𝑓) and 𝛽^(𝑓) can be obtained by replacing (𝛼,𝛽) with (𝛼,𝛽) in (18) and then performing an unconstrained maximization of 𝑙(𝛼,𝛽); for this paper, a standard unconstrained minimization routines available in Matlab software (i.e., fminsearch) was used, setting initial values to (𝛼,𝛽)=(0,0).

Notice that if 𝑍^{(𝑖𝑗)(1)}_SemiP,𝑍^{(𝑖𝑗)(2)}_SemiP,…,𝑍^{(𝑖𝑗)(𝑁)}_SemiP are placed in ascending order and denoted by 𝑍^(𝑖𝑗)_SemiP(1),𝑍^(𝑖𝑗)_SemiP(2),…,𝑍^(𝑖𝑗)_SemiP(𝑁), such that 𝑍^(𝑖𝑗)_SemiP(1)≤𝑍^(𝑖𝑗)_SemiP(2)≤…≤𝑍^(𝑖𝑗)_SemiP(𝑁), then 𝑍^(𝑖𝑗)_SemiP=𝑍^(𝑖𝑗)_SemiP(1) according to (34)—the lowest order statistics (see, for instance, [28]). This fact will be used in estimating the asymptotic behavior of the overall quasi-global semiparametric algorithm, shown in the Appendix.

Notice also that if 𝐖₁ is significantly different from all {𝐖₀^(𝑓)}^𝑁_𝑓=1, then all of the corresponding results {𝑍^{(𝑖𝑗)(𝑓)}_SemiP}^𝑁_𝑓=1 in (35) would yield high values; this outcome would guarantee the lowest order statistics 𝑍^(𝑖𝑗)_SemiP in (34) to hold a high value. Otherwise, if 𝐖₁ is significantly similar to at least one of the samples in {𝐖₂^(𝑓)}^𝑁_𝑓=1, then at least one of the corresponding results in {Z^{(𝑖𝑗)(𝑓)}_SemiP}^𝑁_𝑓=1 would yield a low value; this low value would be assigned to 𝑍^(𝑖𝑗)_SemiP, according to (34).

Since it is unknown a priori whether target spectra are present in 𝐗, the entire 𝐗 needs to be tested. In order to achieve it, all {𝑍^(𝑖𝑗)_SemiP}^{𝑅−𝑛−1,𝐶−𝑛−1}_{𝑖=1,𝑗=1} must be computed according to (34), producing a 2-dim output surface 𝐙^(𝑔)_SemiP. The index 𝑔(1≤𝑔≤𝑀) has been introduced to results produced by (34) in order to denote the repetition number discussed in Section 3.1, yielding 𝐙^(𝑔)_SemiP=𝑍^(11)(𝑔)_SemiP,…,𝑍^{[1(𝐶−𝑛−1)](𝑔)}_SemiP𝑍⋮⋮⋮_{[(𝑅−𝑛−1)1](𝑔)SemiP},…,𝑍^{[(𝑅−𝑛−1)(𝐶−𝑛−1)](𝑔)}_SemiP. (39)

The computation leading to (39) will be independently repeated 𝑀 number of times in order to significantly reduce the probability of contamination (i.e., samples of targets labeled as clutter background). Applying a cutoff threshold to all pixel values 𝑍^{(𝑖𝑗)(𝑔)}_SemiP in 𝐙^(𝑔)_SemiP, such that, pixel values that are above or equal to the threshold become 1 and values below become 0, yielding a binary surface (a probabilistic cutoff threshold for this algorithm is presented in Section 4.3). The 𝑀 binary surfaces are fused using the logic OR operator ⊕, leading to the algorithm’s final output surface 𝐙_SemiP, or𝐙_SemiP=𝑍^(11)(1)_SemiP⊕⋯⊕𝑍^(11)(𝑀)_SemiP,…,𝑍^{[1(𝐶−𝑛−1)](1)}_SemiP⊕⋯⊕𝑍^{[1(𝐶−𝑛−1)](𝑀)}_SemiP𝑍⋮⋮⋮_{[(𝑅−𝑛−1)1](1)SemiP}⊕⋯⊕𝑍^{[(𝑅−𝑛−1)1](𝑀)}_SemiP,…,𝑍^{[(𝑅−𝑛−1)(𝐶−𝑛−1)](1)}_SemiP⊕⋯⊕𝑍^{[(𝑅−𝑛−1)(𝐶−𝑛−1)](𝑀)}_SemiP. (40)

Figure 8 illustrates 𝐙^(𝑔)_SemiP (39) and 𝐙_SemiP (40) through a repeated random sampling and result fusion diagram. The diagram shows M independent paths, where, in each path, independent blocks of data are randomly selected from the input HS data cube so that the entire data cube can be tested against these blocks of data, using a testing window of the same block size. Each path, which is indexed by 𝑔(1≤𝑔≤𝑀), produces a 2-dim output surface (𝐙^(𝑔)_SemiP), where, at the backend of the diagram, all {𝐙^(𝑔)_SemiP}^𝑀_𝑔=1 passes through a logical OR operator on a pixelwise fashion (i.e., only the pixel values at the same pixel location are logically OR’ed), producing a final 2-dim surface 𝐙_SemiP, as shown in (40).

Figure 8 

Quasiglobal, semiparametric anomaly detection algorithm.

The motivation and functionality shown in Figure 8 are summarized as follows: for a given repetition 𝑔(1≤𝑔≤𝑀), assume that the realization of 𝐖₁ from a window location 𝑖𝑗 in 𝐗 is a spectral sample of a target, and the realizations of {𝐖₀^(𝑓)}^𝑁_𝑓=1 are samples of various materials composing the clutter background in 𝐗, that is, the randomly selected blocks of data are not contaminated with target spectra. The semiparametric order statistics in (34) is expected to yield a high value at that 𝑖𝑗 location. Moreover, if the target scale in 𝐗 is larger than 𝑛×𝑛, then the target will be represented by multiple pixels in 𝐙^(𝑔)_SemiP (39), having high values. These pixels are expected to be clustered, hence, accentuating the target spatial location in 𝐙^(𝑔)_SemiP. However, as discussed in Section 3.1, the contamination probability 𝑃(𝑚≥1), for a given 𝑔, increases as a function of increasing 𝑁, see Figure 7. Figure 7 shows further that for a fixed 𝑞, 𝑁 and an adequately large 𝑀, if (for instance) results 𝑍^(22)(1)_SemiP,𝑍^(22)(2)_SemiP,…,𝑍^(22)(𝑀)_SemiP correspond to the same portion of the target at a testing window location (𝑖=2,𝑗=2), then (40) give us the confidence that at least one term in 𝑍^(22)(1)_SemiP,𝑍^(22)(2)_SemiP,…,𝑍^(22)(M)_Semip will have a high value with high probability [1.0−𝑃(𝑚=𝑀)]; after application of a cutoff threshold to results in (39), the high value(s) in reference would be captured by the logic OR operator, for example, (Z^(22)(1)_SemiP⊕⋯⊕Z^(22)(M)_SemiP), as shown in (40) for all 𝑖𝑗 locations. Notice that a target may also be represented by multiple (clustered) pixel locations in 𝐙_SemiP (40).

3.3. Setting the Cutoff Threshold and Other Parameters

For autonomous remote sensing applications, properly setting the algorithm’s parameters is a critical step. This subsection presents a guideline to address this step. For the quasi-global, semiparametric anomaly detection algorithm, the parameters of main concern are 𝑇₁ (the probabilistic cutoff threshold), 𝑁 (the number of randomly selected blocks of data), 𝑀 (the number of testing repetitions), and 𝑞 (the upper bound ratio of target pixels in the data cube over the spatial area of this cube).

Using the asymptotic behavior shown in (A.11) in the Appendix, a cutoff threshold is attained as follow: 𝑇₁=𝑇(𝑎)=𝜇_𝑔+𝑎𝜎_𝑔, (41) where 𝑎=𝑔⁻¹(1−𝜀₁) is the 1−𝜀₁ quantile of 𝑔(𝑧)=𝐿𝑔(𝑧)[1−𝐺(𝑧)]^𝐿−1 (see Appendix), 𝜇_𝑔=_𝑢𝑛₌₁𝑧_𝑢𝑔𝑧_𝑢,𝜎_𝑔=_𝑢𝑛₌₁𝑧²_𝑢𝑔𝑧_𝑢−𝜇²_𝑔 (42) are estimates of the mean and standard deviation, respectively, of the known distribution 𝑔(𝑧)—these estimates can be attained a priori through a simulation experiment using 𝑛 samples of 𝑔(𝑧), and 0≤𝜀₁≤1.

For setting parameters 𝑁 and 𝑀, as discussed in Section 3.1, the quasi-global semiparametric algorithm—ideally—requires an adequately large, which undesirably increases the contamination probability 𝑃(𝑚≥1) per repetition, and an adequately large 𝑀, which desirably decreases the overall cumulative contamination probability 𝑃(𝑚=𝑀). From (24), (25), and (26) and using the log of base 10, a direct transformation leads to [(𝑁=log1−𝑃𝑚≥1)],𝑃(log1−𝑞M=log𝑚=𝑀)log1−1−𝑞^𝑁. (43)(44)

For a given 𝑞, we can fix the values of 𝑃(𝑚≥1), 𝑃(𝑚=𝑀) and obtainand 𝑀 directly using (43) and (44), respectively. As a guideline, 𝑃(𝑚≥1) should be set high, but less than 1.0, so that 𝑁 can also be relatively high and 𝑃(𝑚=𝑀)<1.0; 𝑃(𝑚=𝑀) should be set very low, near zero. As long as the guideline is followed, interestingly, the actual values of 𝑃(𝑚≥1) and 𝑃(𝑚=𝑀) are unimportant. For example, we could fix 𝑃(𝑚≥1)=0.90 and 𝑃(𝑚=𝑀)=0.01, and for 𝑞=0.10, we obtain directly from (43) and (44) parameter values 𝑁≈45 and 𝑀≈44 (Since 𝑁 and 𝑀 are defined as integers, these numbers are rounded off ≈). For the results shown in Section 4, we fixed at once 𝑞=0.10, 𝑃(𝑚≥1)=0.90, and 𝑃_𝑅(𝑚_𝑅=𝑀)=0.015, which by using (43) and (44) yield 𝑁≈22 and 𝑀≈40.

4. Results

The QG-SemiP algorithm was applied to the HS imagery shown in Figure 1, that is, Cube 1 (nadir looking imagery) and Cube 2 through Cube 4 (forward looking imagery), to test for scene (spectral) anomalies. This subsection presents performance summary using results and guideline discussed in Section 3.3 to set algorithm parameters (𝑞,𝑇₁,𝑁,𝑀). Results using forward looking imagery will be discussed first.

4.1. Forward Looking Imagery

Results testing Cube 2 are shown in Figure 9. For display purposes, the output surface 𝐙_SemiP (Figure 9, center and right hand side surfaces) is not shown as a binary surface; instead, each 𝐙^(𝑔)_SemiP(𝑔=1,…,𝑀) output surface is mapped using a pseudocolor map, such that, the brightest pixel values in those surfaces (white colored pixels, representing strongest evidence of anomalies) show the locations of results above or equal to the cutoff threshold 𝑇₁; while other colors (yellow, red, brown, and black) show lesser evidence of anomalies at the corresponding pixel locations, relative to randomly selected blocks of data. (The color black represents no evidence of anomalies.) All of the 𝑀 surfaces 𝐙^(𝑔)_SemiP, using the threshold based colormap in reference, are then summed to yield the output surface shown in Figure 9 through Figure 11. Additional details follow.

(a)

(b)

(c)

Figure 9 

QG-SemiP results testing Cube 2 (a) for scene anomalies; output surface (b) was produced using (𝑞=0.10;𝑇1=𝑇(20);𝑁=3;𝑀=3); output surface (c) was produced using (𝑞=0.10;𝑇1=𝑇(20);𝑁=22;𝑀=40). Bright pixel values (white) in the output surfaces correspond to values above the probabilistic cutoff threshold 𝑇1—depicting the highest confidence level of anomaly presence in the imagery, relative to 𝑁 randomly selected blocks of data. Testing procedure was independently repeated 𝑀 times, as highlighted in Figure 8. Using the available ground truth information of the scene, the white clusters in the far right figure cover about 90% of the motor vehicles (the targets) and no false alarms.

Figures 9(b) and 9(c) show two different outcomes of the quasi-global, semiparametric anomaly detection algorithm, where 𝑛×𝑛 was fixed at once to 20×20 (for all data blocks and testing window sizes) and algorithm’s parameters were set to (𝑞=0.10;𝑇₁=𝑇(20);𝑁=3;𝑀=3)—center display—and (𝑞=0.10;𝑇₁=𝑇(20);𝑁=22;𝑀=40)—right hand side display. The center output surface depicts an example when 𝑁 is not set sufficiently high, hence, obtaining an inadequate representation of the clutter background. In this case, three blocks of data were randomly selected from the scene (most likely from the open field area, since it is the largest area in the scene), and used by the QG-SemiP detector to suppress (according to 𝐙^(𝑔)_SemiP(𝑔=1,…,𝑀) in (39)) the open field in Cube 2, not only once, but most likely 𝑀=3 times. As a result, the three motor vehicles and the canopy area on the upper portion of that scene were accentuated relative to the open field. For this initial experiment, we ignored the binomial distribution model and set parameters 𝑁 and 𝑀 intentionally low and tested Cube 2 to show the undesired result in Figure 9(b).

For parameters accordingly set to (𝑞=0.10;𝑇₁=𝑇(20);𝑁=22;𝑀=40) yielded a significantly better result by detecting only the three motor vehicles in the scene, while suppressing the unknown clutter background, see Figure 9(c). Using the available ground truth information of the scene, the white clusters cover about 90 percent of the motor vehicles’ spatial area (the targets) and no false alarms. As discussed earlier, for many remote sensing applications, targets (if present in the scene) will not cover more than 10 percent of the imagery spatial area. For instance, the motor vehicle shown at broadside in Cube 2 has 25,000 pixels, covering 6.1% of the imagery area (25000/409600). Note that 𝑞=0.1 is robust, because it is independent of targets’ aspect or depression angles, relative to the sensor; independent of the number of targets in the scene; and independent of targets’ scales, relative to other objects in the scene.

The output surface shown in Figure 9(c) shows three manmade objects (motor vehicles) clearly accentuated (pixel values above 𝑇₁) relative to the unknown cluttered environment. It is an achievement, given that no prior information is used about the materials composing the clutter background, or about whether targets are present in the scene, or about targets’ scales relative to other objects in the imagery. But notice in Figure 9 that the standing person in the scene center is not detected, possibly because the window size might be too large and/or there must have some materials in that background (randomly selected) spectrally similar to the materials representing that person (e.g., pants, shirt, skin). Figures 10 and 11 show additional results.

(a)

(b)

Figure 10 

Results testing Cube 3 (a) and corresponding output surface (b). Parameters were set to (𝑞=0.10;𝑇1=𝑇(20);𝑁=22;𝑀=40). A motor vehicle is parked in tree shades—around center of the scene. Using the available ground truth information of the scene, the white clusters cover about 73% of the motor vehicle’s spatial area and no false alarms.

(a)

(b)

Figure 11 

Results testing Cube 4 (a) and corresponding output surface (b). Parameters were set to (𝑞=0.10;𝑇1=𝑇(20);𝑁=22;𝑀=40). A motor vehicle is shown on the left hand side, between top and bottom, behind tree trunks—a sport car. Cube 4 exemplifies a hard case of autonomous clutter suppression. Using the available ground truth information of the scene, the white clusters cover about 44% of the motor vehicle’s spatial area and less than 2% of false alarms.

Figures 10 and 11 show results for Cube 3 and Cube 4, respective; both HS cubes particularly represent difficult cases of clutter suppression. Parameters were also set to (𝑞=0.10;𝑇₁=𝑇(20);𝑁=22;𝑀=40) for both HS cubes.

The guideline described in Section 3.3 for setting parameters also worked very well for both complex scenes depicted in Figures 10 and 11. Both output surfaces clearly accentuates the presence of a motor vehicle—one in tree shades (center in Cube 3) and another motor vehicle parked behind a heavily cluttered environment (center left hand side in Cube 4); see white pixels, or pixel values greater than or equal to 𝑇₁, in both output surfaces in Figures 10 and 11. Setting 𝑇₁=𝑇(20), in essence, means setting the cutoff threshold at 20 standard deviations above the mean of distribution 𝑔(𝑧) in (A.10).

Using the available ground truth information of the scenes in Figure 9 through Figure 11, quantitative comparative performances were obtained via receiver’s operating characteristic (ROC) curves (vertical axis show Pd for probability of detection, and horizontal axis shows Pfa for probability of false alarms) for some of the anomaly detectors mentioned in Section 1. The data cubes were processed using the global RX, k means, GMM, and QG-SemiP. In essence we used the k mean and GMM in place of SemiP in the context of the parallel random sample, but had to take into consideration some of the inherent constraints of these methods. For instance, for k means and GMM, we recorded ROC curves for parameter 𝑁 set to 𝑁=3,5,10,20,and50 and repetition parameter 𝑀=50, while for QG-SemiP 𝑁 was set to 20, 50, 75, and 100 with 𝑀=50. Global RX estimates the mean and covariance from the entire data cube, as discussed in Section 1. Figure 12 shows performance of these detectors using ABC to label QG-SemiP, global to label global RX, KM to label k means, and GMM.

Figure 12 

ROC curve performance testing QG-SemiP (ABC), global RX (Global), k means (KM), and GMM detectors on scene shown in Figure 9.

Although the parameter values for k means and GMM start at a lower value than QG-SemiP, Figure 12 shows that at lower 𝑁 the detectors k means and GMM actually perform much better than for 𝑁 set at higher values, as 𝑁may be interpreted by these algorithms as the number of distinct classes in the scene. Such performance degradation occurs with large 𝑁 because the spatial area that corresponds to each individual 𝑁 block of data is now smaller and samples of the targets, required by the algorithm, need to be included into one of the classes. This outcome contaminated the distribution of the background clutter forcing it to be closer to the distribution of the targets, resulting in performance degradation. The global RX performed reasonably well, as expected since the scene in Figure 9 is relatively less complex. The QG-SemiP detector, also as expected, improved performance as 𝑁 increased.

Performance of these algorithms on the scenes in Figures 10 and 11 are shown in Figures 13 and 14, respectively. Performance degradation of k means, GMM, and global RX are evident from Figure 13, since the target is on tree shades. QG-SemiP performs well for 𝑁>20. In Figure 14, the performance of the k means performed poorly at 𝑁=20, since the target is partially blocked by tree trunks, while GMM performed poorly at 𝑁=3 and 𝑁=5 and was unable to converge at higher 𝑁 values. The global RX surprisingly worked reasonably well, but completely underperforming the QG-SemiP detector.

Figure 13 

ROC curve performance testing QG-SemiP (ABC), global RX (global), k means (KM), and GMM detectors on scene shown in Figure 10.

Figure 14 

ROC curve performance testing QG-SemiP (ABC), global RX (global), k means (KM), and GMM detectors on scene shown in Figure 11.

4.2. Nadir Looking Imagery

For the nadir looking imagery, Cube 1 in Figure 1 (top), ROC curves are also used as a means to quantitatively compare five anomaly detection approaches: local RX, local KRX, local FLD, local SemiP, and QG-SemiP; see Section 1.

Local anomaly detectors process small (𝑛×𝑛) windows of the HS data cube 𝐗, where all the 𝐱_𝑟𝑐(𝑟=1,…,𝑅;𝑐=1,…,𝐶) in 𝐗 are used; modeling is only done at the level of the 𝑛×𝑛 windows, where 𝑛≪𝑅 and 𝑛≪𝐶 (≪ denoting many orders of magnitude smaller than); at the level of the pixel area surrounding these windows. Blocks of data (𝑛×𝑛 windows) that are spectrally different from pixels surrounding them score high using an effective detector in contrast to blocks of data that are not spectrally different from their surrounding pixels. After the detector scores the entire 𝐗, it yields a 2-dim surface 𝐙[𝑎(𝑅−𝑛−1)×(𝐶−𝑛−1) array of scalars], where a cutoff threshold is then compared to the pixel values in 𝐙. Pixels having values greater than the threshold are labeled local anomalies (notice that the SemiP detector will be used in both local and quasi-global versions).

As described in Section 2.1, the set of 14 ground vehicles near the treeline in Cube 1 (Figure 1) constitutes the target set, but, since anomaly detectors are not designed to detect a particular target set, the meaning of false alarms is not absolutely clear in this context. For instance, a genuine local anomaly not belonging to the target set would be incorrectly labeled as a false alarm. Nevertheless, it does add some value to our analysis to compare detections of targets versus nontargets among the different algorithmic approaches.

Figure 15 shows the ROC curves produced by the output of the five algorithms testing Cube 1 for local or scene anomalies. Detection performance was measured using the ground truth information for the HYDICE imagery. We used the coordinates of all the rectangular target regions and their shadows to represent the ground truth target set. As it can be readily assessed from Figure 15, the local SemiP and QG-SemiP anomaly detection approaches outperform the other three techniques on the tested scene, followed by KRX’s performance. The significant display of performance shown in Figure 15 by the semiparamtric algorithms can be further appreciated by taking a closer look at the output surfaces produced by all five detectors, as they show evidences of candidate local and scene anomalies. The intensity of local peaks shown in Figure 16 reflects the strength of the detectors’ evidences. Figure 16 shows that the surfaces produced by FLD, RX, and KRX detectors are expected to be significantly more sensitive (producing false alarms) to changing cutoff thresholds then the ones produced by the local SemiP and QG-SemiP approaches.

Figure 15 

ROC curves for the nadir looking imagery (Cube 1) data scene shown in Figure 1. The semiparametric method, being used in both local and quasiglobal versions, are noticeably less sensitive to different cutoff thresholds; their performances almost achieve an ideal ROC curve for that scene, that is, a step function starting at point (PFA=0,PD=1).

Figure 16 

Detection algorithms’ output surfaces for Cube 1 (far left). The intensity of local peaks reflects the strength of evidences as seen by different anomaly detection approaches. Boundary issues were ignored in this test; surfaces were magnified to about the size of the original image only for the purpose of visual comparison. FLD, RX, KRX, and SemiP performed local anomaly detection by testing spectra within a testing window (red square shown in the scene display—far left, top) to spectra surrounding the testing window (outer window bounded by yellow lines). QG-SemiP performed global anomaly detection, as presented in this paper.

Spatial areas shown in Cube 1 containing the presence of clutter mixtures (e.g., edge of terrain, edge of tree clusters), where FLD, RX, and KRX yield a high number of potential false alarms (false anomalies), are suppressed by the SemiP approach, local, and quasi-global. The reason for this suppression is that, as part of the overall comparison strategy, the semiparametric model combines both the test and reference samples in order to estimate the underlying PDF of the reference sample, as shown by simulation earlier. As such, the semiparametric test statistic ensures that a component of a mixture (e.g., shadow) will not be detected as a local anomaly when it is tested against the mixture itself (e.g., trees and shadows). Performances of such cases are represented in Figure 16 in the form of softer anomalies (significantly less-dominant peaks) in the local SemiP’s and QG-SemiP’s output surfaces. It is evident from Figure 16 that both versions of the SemiP detectors perform remarkably well accentuating the presence of dominant local or scene anomalies (e.g., targets) from softer anomalies (e.g., transitions of distinct regions). The natural ability of the semiparametric model to manage spectral mixture can best explain the local SemiP and QG-SemiP superior ROC-curve performances shown in Figure 15.

5. Conclusion

This paper introduced an adaptive quasi-global, semiparametric anomaly detection algorithm and evaluated the approach using real HS imagery, where targets (manmade objects) are found in difficult natural clutter backgrounds viewed from two different perspectives—nadir and forward looking. The algorithm features a semiparametric test statistic, which has been recently found to be robust against spectral mixture, and applies random sampling of the imagery to test for anomalies. Random sampling and testing are repeated a number of times in order to mitigate the probability of contamination (spectral samples of candidate targets being sampled and used as clutter reference samples). As such, the algorithm requires no prior information (e.g., a spectral library of the clutter background and/or target, target size, or shape). The algorithm is free from training requirements. We found that the semiparametric model has a natural ability to handle mixtures, although an exhaustive survey of the literature reveals that this fact has never been noticed before by practitioners of the model in other fields of study (e.g., biotechnology).

The repeated sampling and testing procedure was modeled by the binomial family of distributions, where the only target related parameter 𝑞 (the upper bound proportion of target pixels potentially covering the spatial area of the imagery) is robust—thus invariant—to different sizes and shapes of targets, number of targets present in the scene, target aspect angle, partially obscured targets, or sensor viewing perspective. The paper also discussed how other parameters (𝑁, the total number of sampled data blocks to take from the HS imagery, and 𝑀, the number of process repetitions) can be automatically set using a simple guideline.

The algorithm fuses intermediate results through the application of minimum order statistics and logic OR operation. The paper presented the algorithm’s asymptotic behavior under the null hypothesis, when either the null or the alternative hypothesis is true, for the two-sample test case and the multi-sample test case, where the cumulative probability of the algorithm making mistakes was derived. Using the cumulative probability, a cutoff threshold can be determined from a user specified probability of error. This is a desired feature giving the user some control of predicted errors.

The inherent challenges in adequately modeling spectral variability of targets, while managing spectral variability between targets, have prompted the introduction of a more robust family of algorithms that attempts to detect targets as being anomalous to an unknown natural clutter background. This paper presented an approach that inherently addresses many of the problems and issues pertaining to anomaly detection applications. The advantages of using anomaly detection algorithms have been discussed in this paper; however, it should be emphasized that targets would not be detected as specific manmade objects; they would be detected as anomalies. This is a limitation.

Appendix

Asymptotic Behavior of the Quasiglobal Semiparametric Algorithm

This section shows an analytical asymptotic analysis of the quasi-global semiparametric anomaly detection algorithm. In particular, we would like to investigate the algorithm’s cumulative probability of rejecting the null hypothesis, when either the null or the alternative hypothesis is true; that is, the algorithm’s probability of making mistakes. We will look first at the asymptotic behavior of the two sample test case, where the detector tests a reference sample against a test sample; then we will look at the multisample test case, which uses order statistic to reduce 𝑁 results to a single result.

Two-Sample Test (2ST)
In order to declare an anomaly, a decision threshold 𝑇 must be chosen; hopefully, separating without errors the null and the alternative hypotheses in some decision space. In the paper’s context, values of 𝑍_SemiP in (23) greater than 𝑇 are automatically labeled as anomalies. And since, in real world applications, decision errors are unavoidable, we would like to know whether the asymptotic behaviors of these errors can be determined, and whether they are favorable. The power function can provide those answers for the two-sample test (2ST) case. Using (13), the power function of the semiparmetric test statistic for the two-sample test (2ST) case is as follows: 𝜓𝛽=𝑃_𝛽=0𝑍_SemiP𝑃>𝑇_𝛽≠0𝑍_SemiP>𝑇. (A.1)

In essence, the power function 𝜓 yields the cumulative probability 𝑃 of rejecting the null hypothesis 𝐻₀ in (13) when either 𝐻₀ (𝛽=0) or the alternative 𝐻₁ (𝛽≠0) is true. This rejection region is 𝑍_SemiP>𝑇, where 𝑍_SemiP is defined in (23) and 𝑇 is the decision threshold. Notice in (A.1) that 𝜓 under 𝐻₀ corresponds to the well known type I error, or the probability of missing (i.e., the probability of rejecting 𝐻₀, given that 𝐻₀ is true) and that 𝜓 under 𝐻₁ corresponds to the complement of the type II error, or probability of false alarms (i.e., one minus the probability of rejecting 𝐻₁, given that 𝐻₁ is true). The type I and type II errors constitute the only error types encountered in the context of our discussion. In the ideal case, 𝜓 yields 0.0 when 𝐻₀ is true and 1.0 when 𝐻₁ is true. Except in trivial situations, this ideal cannot be attained. So, one of our goals is to show that 𝜓 tends in probability to 𝜀 (a scalar controlled by the user), when 𝐻₀ is true, and that 𝜓 tends in probability to 1.0, when the alternative hypothesis 𝐻₁ is true.

If 𝐻₀ in (13) is true, the semiparametric detector has the asymptotic behavior shown in (23), and the type-I error is readily obtained by the following: 𝑃_𝛽=0𝑍_SemiP>𝑇→−_𝑛→∞𝑃(𝜉>𝑇)=𝜀, (A.2) where 𝜉 is a Chi square distributed random variable with 1 degree of freedom, 𝑍_SemiP as defined in (23), 𝜀 and 𝑇 are nonnegative real values.

Setting 𝜓(𝛽)=𝑃_𝛽=0(𝑍_SemiP>𝑇), 𝜓 is indeed an asymptotic size 𝜀 test, which is controlled by the user.

Now consider an alternative parameter value, such that 𝛽≠0, and let 𝜁² be the true variance in (21) and 𝜁̂² be a estimator of 𝜁², or 𝜁²=𝜌⁻¹1+𝜌²𝑣²,𝜁̂²=𝜌⁻¹1+𝜌²𝑣², (A.3) where𝜌=𝑛₁/𝑛₀ and 𝑣² is defined in (22).

From (23), we can now write 𝑍_SemiP=𝛽−𝛽𝜁²/𝑛_𝐴+𝛽𝜁²/𝑛_𝐵²𝜁²𝜁²_𝐶̂. (A.4)

Notice in (A.4) that, as 𝑛₁ and 𝑛₀ go to +∞, 𝜌=𝑛₁/𝑛₀ tends to 1 and 𝜁² tends to 𝜁² for 𝑣²>0. According to (21), the term 𝐴 in (A.4) converges in distribution to the standard Normal, 𝑁(0,1), as 𝑛₁ and 𝑛₀ (hence, 𝑛) go to +∞, no matter what the values of 𝛽 or 𝜁² are. Note also that the term 𝐵 converges to +∞ or −∞ in probability, as 𝑛 goes to +∞, depending on whether 𝛽 is positive or negative. Since the estimator 𝑔₀ in (38) has been shown [21] to be biased, as 𝑛₁ and 𝑛₀ go to +∞, 𝜁̂² tends to a constant, see definition of 𝑣² in (22); leading the term 𝐶 also to a constant, no matter what values of 𝜁² is. Thus, 𝑍_SemiP converges to +∞ in probability and 𝑃_𝛽≠0reject𝐻₀=𝑃_𝛽≠0𝑍_SemiP>𝑇→−_𝑛→∞1. (A.5) In this way, the semiparametric test statistic shown in (23) also has the properties of asymptotic size 𝜀 and asymptotic power 1, which is highly desired.

Multisample Test (mST)
The discussions in Sections 2.3 and 4.1 ensure that the output of the semiparametric 2ST has two asymptotic outcomes: Z_SemiP→−_𝑛→∞𝜒²₁ in distribution, if 𝐻₀ in (13) is true, or 𝑍_SemiP→−_𝑛→∞+∞ in probability, if 𝐻₁ is true.

Using results leading to (13) and the order statistics 𝑍^(𝑖𝑗)_SemiP=min_{1≤𝑓≤𝑁}𝑍^{(𝑖𝑗)(𝑓)}_SemiP in (34), for the multi-sample test (mST) case, we propose the following null 𝐻₂ and alternative 𝐻₃ hypotheses 𝐻₂∶atleastone𝛽^(𝑓)𝑁_𝑓=1𝐻=0,₃∶all𝛽^(𝑓)𝑁_𝑓=1≠0, (A.6) where 𝛽^(𝑓) is the true logistic function parameter—see (12)—corresponding to the 𝑓th randomly selected block of data from a HS data cube 𝐗.

Now, consider the following: for a given spatial location in 𝐗, let 𝑍^(𝑓)_SemiP be the semiparametric detector’s output for the 𝑓th block of data, and assume, without loss of generality, that each one of the first L outputs in the independent sequence of results (1≤𝐿≤𝑁, where 𝑁 is the total number of randomly selected blocks of data in 𝐗) has the asymptotic chi-square behavior shown in (23), and that each one of the remainder results has the asymptotic behavior tending to +∞, or 𝑍_SemiP𝐻2𝑍=min^{(1)SemiP𝑛→∞}−−−−−→𝜒²₁𝑍−−−−−⁽²⁾_SemiP^𝑛→∞→𝜒²₁⋮𝑍−−−−−^(𝐿)_SemiP^𝑛→∞→𝜒²₁𝑍^(𝐿+1)_SemiP^𝑛→∞⋮𝑍−−−−−→+∞_{(𝑁)SemiP𝑛→∞}−−−−−→+∞. (A.7) Under the null hypothesis 𝐻₂ in (A.6), 𝑍_SemiP in (A.7) is bounded because, as 𝑛→∞, 𝑍_SemiP will converge in law to the distribution of the lowest order statistics. (The order statistics of a random sample 𝑍₁,…,𝑍_𝑁 are the sample values placed in ascending order. They are often denoted by 𝑍₍₁₎,…,𝑍_(𝑁), where 𝑍₍₁₎=min_{𝑓≤𝑓≤𝑁}𝑋_𝑓 and 𝑍_(𝑁)=max_{𝑓≤𝑓≤𝑁}𝑋_𝑓.) To attain an approximation of the type I error using (A.7), we first ignore all the components in (A.7) that converge in probability to +∞, then we consider only the components that converge in distribution, that is, (𝑍⁽¹⁾_SemiP,𝑍⁽²⁾_SemiP,…,𝑍^(𝐿)_SemiP). The distribution of 𝑍_SemiP(1)=min_{𝑓≤𝑓≤𝐿}𝑍^(𝑓)_SemiP (A.8) from the culled sequence can be attained with the application of Theorem 1.

Theorem 1. Let 𝑋₍₁₎,…,𝑋_(𝑛) denote the order statistics of a random sample from a continuous population with cumulative distribution function (cdf) 𝐹(𝑥) and pdf 𝑓(𝑥). Then the pdf of 𝑋_(𝑗) is 𝑓(𝑥)=𝑛!!!𝑓(𝑥]𝑗−1𝑛−𝑗)[𝐹(𝑥)^𝑗−1[1−𝐹(𝑥)]^𝑛−𝑗, (A.9) where (⋅)! denotes the factorial operator.

The proof of Theorem 1 can be found in [28].

Using Theorem 1 with 𝑗=1 and 𝑛=𝐿, the pdf of 𝑍_SemiP=𝑍_SemiP(1) under 𝐻₂ in (A.6) is 𝑔(𝑧)=𝐿𝑔(𝑧)[1−𝐺(𝑧)]^𝐿−1, (A.10) where 𝑔(𝑧) is the Chi square pdf with 1 degree of freedom and 𝐺(𝑧) is the corresponding cdf.

Denote the 𝑘th logistic function parameter 𝛽^(𝑘) in (A.6) to correspond to the one of the minimum order statistics 𝑍_SemiP(1). As the sample size increases in 𝑍_SemiP(1), that is, 𝑛=𝑛^(𝑘)→∞, the probability of rejecting the null hypothesis 𝐻₂ in (A.6), when 𝛽^(𝑘)=0, converges to 𝜓𝛽^(𝑘)=𝑃_𝛽^(𝑘)₌₀𝑍_SemiP>𝑇₁→−_𝑛→∞𝑃𝜉>𝑇₁=𝜀₁, (A.11) where 𝜉 is a random variable distributed by 𝑔(𝑧), as defined in (A.10); 𝑇₁ a nonnegative real value; 𝜀₁ is a positive real value, controlled by the user.

The variable 𝜓 in (A.11) is the type I error under 𝐻₂ for the mST case, and it is indeed an asymptotically size 𝜀₁ test.

Now consider the alternative hypothesis 𝐻₃ in (A.6), where all {𝛽^(𝑓)}^𝑁_𝑓=1≠0. From (A.7) one can write 𝑍_SemiP𝐻3𝑍=min⋮^{(1)SemiP𝑛→∞}𝑍−−−−−→+∞_{(𝑁)SemiP𝑛→∞}−−−−−→+∞. (A.12)

From (A.12), 𝑍_SemiP will converge in probability to +∞, hence, the probability 𝑃 of rejecting the null hypothesis 𝐻₂, given that 𝐻₃ is true, tends to 1.0, or 𝑃_𝛽^(𝑘)_≠0rejecting𝐻₂=𝑃_𝛽^(𝑘)_≠0𝑍_SemiP(1)>𝑇₁→−_𝑛→∞1. (A.13)

In this way, the quasi-global semiparametric anomaly detection algorithm has the desired properties of asymptotic size 𝜀₁, which is controlled by the user, and asymptotic power 1.0.

Notations

 Bold upper case letters may denote a data cube (3 dimensions) or a matrix (e.g., 𝐗,𝐖₀), where the specific case in use is defined in the text. Lower cases letters denote vectors (bold) or sequences (not bold) (e.g., 𝐱, 𝑥₁=(𝑥₁₁,…,𝑥_1𝑛1)). PDF or pdf: Probability density function. IID or iid: Independent and identically distributed. 𝐑^𝑑1×𝑑2-𝑑₁ by 𝑑₂ dimensional set of real numbers.  ∈ denotes set belonging  HS: Hyperspectral 𝐗: Observed hyperspectral data cube with dimensions of 𝑅 rows, 𝐶 columns, and 𝐾 bands. 𝐱_𝑟𝑐: Observed spectrum contained in 𝐗 with spatial indexes (𝑟=1,…,𝑅) and (𝑐=1,…,𝐶). A slideing 𝑛×𝑛 window is a 3-dim subset of 𝐗, containing 𝑛⋅𝑛 spectra.  𝐖₁∈𝐑^𝐾×𝑛1: A matrix representing a hyperspectral sample being observed from a sliding 𝑛×𝑛 window in 𝐗 (also referred to herein as a test sample); this is a rearranged version of a 3-dim subdata cube, where vertical direction is the dimension 𝐾 of bands and the horizontal direction is the dimension of countable samples with sample size 𝑛₁=𝑛².  𝐲_1ℎ∈𝐑^𝐾(ℎ=1,…,𝑛₁): An observed spectrum of 𝐾 bands contained in 𝐖₁. 𝑔₁(𝐲|𝜽): Multivariate joint PDF of 𝐲₁₁,…,𝐲_1𝑛1. 𝐖₀∈𝐑^𝐾×𝑛0: A matrix representing a hyperspectral sample labled as a reference sample of sample size 𝑛₀, having the same specifications of 𝐖₁ except perhaps the sample size (𝑛₀ may be different from 𝑛₁).  𝐲_0ℎ∈𝐑^𝐾(ℎ=1,…,𝑛₀): An observed spectrum of 𝐾 bands contained in 𝐖₀. 𝑔₀(𝐲|𝜽): Multivariate joint PDF of 𝐲₀₁,…,𝐲_0𝑛0. ∇₀∈𝐑^{(𝐾−1)×𝑛0}: Output from differentiating 𝐖₀. ∇₁∈𝐑^{(𝐾−1)×𝑛1}: Output from differentiating 𝐖₁. 𝑥₀=(𝑥₀₁,𝑥₀₂,…,𝑥_0𝑛0): Univariate sequence used as the reference sample. 𝑥₁=(𝑥₁₁,𝑥₁₂,…,𝑥_1𝑛1): Univariate sequence used as the test sample. 𝑔₀(𝑥): Univariate PDF labeled as reference. 𝑔₁(𝑥): Univariate PDF labeled as test 𝑡=(𝑥₁₁,…,𝑥_1𝑛1,𝑥₀₁,…,𝑥_0𝑛0)≡(𝑡₁,…,𝑡_𝑛): Sample concatenation, combining samples. 𝑔₀(𝑡): estimator of 𝑔₀(𝑥).  𝑔₀(𝑡): estimator of 𝑔₀(𝑡).  𝐻_𝑖: Statistical hypothesis 𝑖. 𝑍_SemiP: Univariate output of the semiparametric detector. 𝑃: Cumulative probability function, using the binomial family of PDFs as base PDF. N: The number of randomly selected blocks of data used to represent background objects. Trial (or process): Take 𝑁 random blocks of data from the data cube under test, label them as reference background objects, and—using the semiparametric detector and a sliding window across the data cube—test the entire data cube against the same set of 𝑁 reference blocks of data. 𝑀: The number of trials (repetitions or parallel processes). 𝑃_𝑔(𝑚≥1): Cumulative probability of contamination, that is, probability of labeling a randomly selected target sample as a background sample at the 𝑔th trial or process. 𝑃: Overall cumulative probability that all of the trials (or processes) are contaminated with at least a contaminated sample from the randomly selected set of reference samples {𝐖₀^(𝑓)}^𝑁_𝑓=1. 𝑍^(𝑖𝑗)_SemiP: Retains the lowest order statistics from a set of 𝑁 semiparametric detector’s results. 𝐙^(𝑔)_SemiP: The 2-dimensional output surface, consisting of 𝑍^(𝑖𝑗)_SemiP values, from the 𝑔th trial. 𝑇₁: Adaptive cutoff threshold for 𝐙^(𝑔)_SemiP. 𝐙_SemiP: A final binary 2-dimensional output surface of the quasi-global semiparametric detector.