Multiscale Superpixel Guided Discriminative Forest for Hyperspectral Anomaly Detection

Abstract

Recently, the isolation forest (IF) methods have received increasing attention for their promising performance in hyperspectral anomaly detection (HAD). However, limited by the ability of exploiting spatial-spectral information, existing IF-based methods suffer from a lot of false alarms and disappointing performance of detecting local anomalies. To overcome the two problems, a multiscale superpixel guided discriminative forest method is proposed for HAD. First, the multiscale superpixel segmentation is employed to generate some homogeneous regions, and it can effectively extract spatial information to guide anomaly detection for the discriminative forest in local areas. Then, a novel discriminative forest (DF) model with the gain split criterion is designed, which enhances the sensitivity of the DF to local anomalies by the utilization of multi-dimension spectral bands for node division; meanwhile, the acceptable range of hyperplane attribute values is introduced to capture any unseen anomaly pixels that are out-of-range in the evaluation stage. Finally, for the high false alarm rate situation in the existing IF-based algorithms, the multiscale fusion with guided filtering is put forward to refine the initial detection results from the DF. In addition, the extensive experimental results on four real hyperspectral datasets demonstrate the effectiveness of the proposed method.

1. Introduction

Hyperspectral images (HSIs) have abundant spectral information of geographical materials, which is different from the traditional optical and spectral images, and this feature makes different ground targets more discriminative. Accordingly, the HSI derives many research directions, such as target detection [1,2], anomaly detection [3,4], band selection [5,6], hyperspectral unmixing [7], hyperspectral denoising [8,9], etc.

Hyperspectral anomaly detection (HAD) [10,11] is a challenging task in hyperspectral image (HSI) processing, and its key purpose is to capture anomaly targets from the background by the utilization of spectral and spatial differences with no prior knowledge. Benefiting from the above characteristic, the HAD is widely used in civil and military applications.

The isolation forest (IF) methods [12] are attracting more attention to the HAD in the last three years, which take advantage of the fact that anomalies are “few and different” in a given HSI, thus anomaly pixels are more susceptible to be isolated than background pixels during the construction processing of an isolation forest. To be specific, Li et al. [13] propose the kernel isolation forest detection (KIFD) method, which combines the kernel principal component analysis (KPCA) with the IF to obtain the initial detection maps, then the iteration strategy is applied to refine the detection results. Afterward, a novel isolation forest method based on multiple features (MFIFD) [14] is proposed for HAD. The multiple features including the Gabor, extend morphological profile (EMP), extend multi-attribute profile (EMAP) and spectral features are applied to build isolation forest to improve the detection accuracy. To better define a pixel as an anomaly or background, a subspace selection-based discriminative forest (SSDF) [15] is put forward.

Although the above IF-based methods have achieved a good effect in the HAD, these methods still have two deficiencies, i.e., the poor performance of detecting local anomalies and a high false alarm rate (FAR), due to the limited utilization capacity of HSI spectral-spatial information. The two problems are analyzed from the angles of methodology and detection results. In the methodology, the IF-based methods usually adopt the path length in the global samples to measure the anomaly degree of test pixels [16], but the essence of anomalies in an HSI is that the anomaly degree of each pixel depends on the difference among its surrounding pixels. By comparison, the global way in the IF-based methods does not correspond with the anomaly essence and ignores the discriminability of anomalous pixels in local spatial areas. Furthermore, using an individual spectral band is not very effective in each node partition of isolation trees. Consequently, these two points lead to poor performance in detecting local anomalies. On the other hand, the IF-based methods in HSI complex scenes ignore the importance of spatial features and cannot distinguish the background pixels well, resulting in suffering from a lot of false alarms. From the perspective of detection results, the anomaly detection maps are shown in Figure 1, and the SSDF and KIFD methods are the isolation forest methods. It is obvious that the SSDF method has poor performance in detecting local anomalies. To be specific, as is shown in Figure 1c, the SSDF method is not able to detect well anomaly targets which are several buildings in the upper left corner. In addition, the SSDF and KIFD methods always erroneously detect the green background pixels as anomalies, causing a lot of false alarms, especially the KIFD method.

In order to solve the two problems, the multiscale superpixel guided discriminative forest (MSGDF) algorithm is proposed in Figure 2. The major contributions are as follows.

(1)

The strategy of multiscale superpixel guided discrimination forest is proposed for hyperspectral anomaly detection in this article. Concretely, the multiscale superpixel segmentation is adopted to mine HSI spatial information and make anomalies more discriminative in local areas, which can better guide the detection of local anomalies for DF.

(2)

A novel discriminative forest (DF) is designed. On the one hand, a gain split criterion based on multi-dimension spectral bands is utilized to make the DF more sensitive to anomaly pixels; On the other hand, the acceptable range of hyperplane attribute values is introduced to capture any unseen anomaly pixels that are out-of-range in the evaluation stage.

(3)

Aiming at the problem of high false alarm rates in the existing IF-based methods, the multiscale fusion with guided filtering (MFGF) module is put forward in the paper, which integrates multiscale detection maps with different spatial features to refine the initial DF results.

(4)

Compared with the state-of-the-art algorithms, the experimental results on four HSI datasets demonstrate the superiority of the proposed method.

The remainder of this paper is organized as follows. The details of the proposed MSGDF method are described in Section 2. Section 3 introduces the experiments and results. The discussion and conclusion are presented in Section 4 and Section 5, respectively.

2. Related Work

Many HAD methods have been developed in the past few decades [10]. One of the most classical algorithms is the Reed–Xiaoli (RX) [17], which proposes the hypothesis that the background features obey the Gaussian multivariate distribution, and the mean vector and covariance matrix of the whole HSI pixels are employed to represent the background, then the Mahalanobis distance is applied to measure the anomaly degree of test pixels. Whereas some scholars believe that the assumption of RX is flawed in the real HSIs, thus several variants are designed to optimize the RX, such as the Kernel RX (KRX) [18], the local RX (LRX) [19], weighted RX (WRX) [20] and the other RX-based methods [21,22].

Subsequently, the reconstruction-based theories have aroused great concern, and their intention is to reconstruct all HSI pixels by the feasible envisaged models. Reconstruction-based algorithms include the low-rank representation (LRR) methods [9,23], the collaborative representation (CR) methods [24,25], and the deep learning methods [26,27,28]. Among LRR methods, the low-rank and sparse representation (LRASR) [9] is the first model to apply the LRR to HAD, and a linear combination of background dictionary’s atoms is employed to denote each background pixel and separates the anomalies. For optimizing the local features and spatial information. The concept of graph regularization and total variation is introduced into the LRR model [23]. Meanwhile, in the aspect of describing different noise characteristics of the LRR, the low-rank and sparse decomposition model with the mixture of the Gaussian model (LSDM-MoG) [29] is presented for the HAD, which adopts the mixture Gaussian to represent the sparse component consisting of noise and anomalies in an HSI. In CR methods, the background pixels can be indicated by their adjacent atoms while the anomaly pixels cannot [24]. Based on this idea, Tu et al., put forward the modified CRD method, which combines the density peak clustering model and CR detector to improve the detection performance.

Additionally, with the widespread application of deep learning [30,31,32], the generative adversarial network (GAN) and autoencoder (AE) methods are employed to construct the background to get the response of anomalies by residual between the HSI and background [26]. For example, the autoencoder adversarial network (AEAN) [28] utilizes the AE of three different dimensions to generate the reconstructed map, then a weighted RX method is to acquire the detection result. Due to suffering from the contamination of anomaly pixels in the processing of background distribution estimation of AE or GAN, a novel semi-supervised GAN method is designed to deal with this problem [27].

3. Methodology

This section is divided into three parts: feature extraction with multiscale superpixel segmentation, discriminative forest, and multiscale fusion with guided filtering. First, the feature extraction with superpixel segmentation is to obtain multiscale local spatial information features by setting different numbers of superpixels, which guides the building of discriminative forest in different blocks of an HSI. Then, based on the above multiscale superpixel features, discriminative forests are built to evaluate the anomaly score of each pixel and generate multiscale initial detection maps. It is worth noting that each superpixel block in an HSI is a discriminative forest to train and test. Finally, the multiscale fusion with guided filtering is to integrate multiscale initial detection maps to obtain the final detection result, which can effectively optimize the detection rate and false alarm rate.

3.1. Feature Extraction with Multiscale Superpixel Segmentation

Based on the fact that anomalous objects occupy very few pixels in an HSI, simple linear iterative clustering (SLIC) [33] is adopted to guide anomaly detection in the local areas by generating some local homogeneous areas containing anomaly pixels. The way fits the definition of HSI anomalies [34,35,36], and it makes anomaly pixels discriminative in local areas. Furthermore, the good HAD effect in the spectral–spatial anomaly detector based on improved isolation forest (SSIIF) [16] confirms the importance of superpixel segmentation for IF methods. Therefore, the SLIC is applied to mine local spatial features in this paper. However, the difference from the SSIIFD method is that we utilize multiscale superpixel maps, which can extract HSI spatial information better than singlescale superpixel maps of the SSIIFD method.

Specifically, a 3D HSI can be first transformed into a 2D matrix $X\in {ℝ}^{N\times D}$, where N means the pixel number and D denotes the number of spectral bands.

Then, the 2D matrix is fed into the principal component analysis (PCA) model [37] to acquire the top three principal components $Y\in {ℝ}^{N\times 3}$. The PCA model can be expressed as

$$\underset{P^{T}P=I}{\min }\hspace{0.17em}tr(P^T{X}^{T}XP)$$

(1)

where $P\in {ℝ}^{D\times 3}$ denotes a projection matrix; $I\in {ℝ}^{3\times 3}$ is an identity matrix; $tr(\cdot )$ indicates the trace of a matrix. The top three principal components $Y\in {ℝ}^{N\times 3}$ are formulated as

$$Y\in {ℝ}^{N\times 3}=(X-E(X))P$$

(2)

where $E(\cdot )$ is the mean function.

Finally, the features of PCA are imported to the SLIC model. By setting different numbers of superpixels, the SLIC can generate multiscale superpixel maps to index multiscale local HSI samples. The SLIC model can be defined as

$$\begin{array}{l}{d}_c=\sqrt{{(l_j-l_i)}^2+{(a_j-a_i)}^2+{(b_j-b_i)}^2}\\ d_s=\sqrt{{(x_j-x_i)}^2+{(y_j-y_i)}^2}\\ D^{\prime }=\sqrt{{(\frac{d_c}{N_c})}^2+{(\frac{d_{\mathrm{s}}}{N_s})}^2}\end{array}$$

(3)

where i and j represent the two different pixels in the features of Y; $d_c$ and $d_s$ are three spectral distances and the spatial distance between the pixel i and j, respectively; $D^{\prime }$ is defined as the final distance; $N_c$ denotes the coefficient of max spectral distance, and it is set to 10 in this paper; and $N_s$ denotes the coefficient of max intra-class spatial distance, which is formulated as

$$N_s=\sqrt{M/K}$$

(4)

where M denotes the number of the input HSI pixels, and K means the number of superpixels setting.

Each superpixel block sample is denoted by ${\mathrm{X}}_w^{\mathrm{r}}$, where r and w are the number of scales and the number of superpixels, respectively. The r-th scale HSI feature ($X^r$) is defined by ${[X_1^r]}_{N_1^r\times D},\hspace{0.17em}{[X_2^r]}_{N_2^r\times D},\hspace{0.17em}\cdots ,\hspace{0.17em}{[X_w^r]}_{N_w^r\times D}$, where $N_1^r{+N}_2^r+\cdot \cdot \cdot {+N}_w^r=N$, i.e., $X_1^r\cup \hspace{0.17em}{X}_2^r\hspace{0.17em}\cup \cdot \cdot \cdot \cup X_w^r=X^r$. As shown in Figure 1, each HSI superpixel block of ${\mathrm{X}}_w^{\mathrm{r}}$ are taken as a subsample to build the local discriminative forest individually, and discriminative forests of the w HSI superpixel block make up an HSI detection map.

3.2. Discriminative Forest

The discrimination forest consists of two phases, namely the training stage and the evaluation stage.

3.2.1. Training Stage

In the training stage, the Diforest is an ensemble of t trees, and the construction of each tree is described in Algorithm 1. Each tree subsample (X^′) is randomly selected from n pixels of an HSI superpixel block ($X_w^r$), and $\hspace{0.17em}{\mathbf{X}}^{\prime }\subset X_w^r$, |X^′| = n. An input pixel is defined as x = {x₁, …, x_D}. Subsequently, the subsamples are recursively divided into two sub-nodes by the best-split threshold, until tree height H reaches the limit H_max ($H_{\max }{=log}_{2}n$) or pixel number of a sub-node is not more than two.

Algorithm 1: Building a tree in Diforest (X′, q, H, Hmax)

Input: X′—input data, q—number of attributes selected from D spectral bands, H—current tree height, which is initialized to zero, Hmax—height limit

Output: an iTree T

1:    if |X′| ≤ 2 or H ≥ Hmax then

2:       return exNode {Size ← |X′|}

3:    else

4:       Z ← generate a new attribute set by q, X′ and Formula (5)

5:       S ← obtain the splitting threshold by Z and Formula (6)

6:       Xl′ ← {z∈Z|zi ≤ S}

7:       Xr′ ← {z∈Z|zi > S}    8:       ν ← {ν = zmax − zmin |z∈Z|}

9:       return inNode{Left ← iTree (Xl′, q, H + 1, Hmax)

10:          Right ← iTree (Xr′, q, H + 1, Hmax)    11:          SplitPlane ← Z    12:          SplitThreshold ← S

13:          Upper Limit ← S + ν

14:          Lower Limit ← S − ν

15:    end if

The IF-based methods of HAD usually employ a single spectral band in every node division, which is not very effective to exploit HSI spectral information [38]. To better utilize spectral features, a novel split threshold is introduced in this paper, and it is obtained by the gain criterion in a hyperplane of q randomly selected spectral bands. First, we construct a split hyperplane f that is non-axis-parallel to the original HSI properties, and each pixel of X^′ is defined as follows:

$$f(x)={\displaystyle \sum _{j\in Q}{c}_j\frac{x_j-\mu (X_j^{\prime })}{\sigma (X_j^{\prime })}}.$$

(5)

where Q has q attribute indices which are randomly selected from D spectral bands in an HSI; c_j is a random coefficient of [–1, 1]; $\mu (\cdot )$ is the average function and $\sigma (\cdot )$ is the standard deviation function; $X_j^{\prime }$ is the jth spectral band attributes from X^′, and x_j is the jth band attribute in an input pixel of x. All pixels of X^′ are transformed into a new attribute set (Z) by Equation (5). Then, based on the Z, the gain splitting criterion S_gain is adopted to obtain the best split threshold S, and it is defined as

$$S_{gain}=\frac{n_{Z_{\mathrm{l}}}{\sigma }_{Z_{\mathrm{l}}}+n_{Z_r}{\sigma }_{Z_r}}{n}\hspace{0.17em}(\hspace{0.17em}{Z}_l=\{z_i\le S\},\hspace{0.17em}{Z}_r=\{z_i>S\}\hspace{0.17em}),$$

(6)

where ${\mathrm{Z}}_{\mathrm{l}}{\displaystyle \cup }{\mathrm{Z}}_{\mathrm{r}}=\mathrm{Z}$, Z_l and Z_r are the samples of a left node and the samples of a right node, respectively; n denotes the total pixel number of Z. In Equation (6), the idea that the sum of within-class standard deviation is minimized to clearly separate two different distributions in each split is utilized. In other words, when the gain splitting criterion S_gain within two sub-nodes is minimized, the best division threshold S is obtained.

3.2.2. Evaluation Stage

The evaluation stage is described in Algorithm 2. Each test pixel x is input into the trained Diforest to assess its anomaly score by the path length, which is the edge number from the root node to the external node in an iTree. If the path length of a test point x is shorter, its anomaly degree is higher, conversely, the anomaly degree is lower. When x is judged to be an external node in an iTree, the value of c(T.Size) is employed to evaluate the path length of the sub-node of unsuccessful tree construction (|X^′| ≤ 2 or h ≥ H_max). The c(T.Size) can be formulated as:

$$c(m)=2H(m-1)-2(m-1)/m,$$

(7)

where H(⋅) is the harmonic number which is represented by ln(i) + 0.5772156649 (Euler’s constant) [39], and m means the sample number of a sub-node. When m is 1 or 2, the value of c(m) is 0; when m is set to the total subsample size n, c(n) denotes the average height of a tree; if m is set to the size of Step 2 in Algorithm 2, c(T.Size) is as the adjust path length.

Algorithm 2: Path Length (x, T, h)

Input: x—a test pixel, T—an iTree, h—number of edges from the root node to the external node, which is initialized to zero

Output: path length of x

1:    if T is an exNode then

2:        return h + c(T.Size)

3:    end if

4:    z ← T. SplitPlane(x)

5:    S ← T.SplitThreshold

6:    if z > S then

7:       return Path Length(x, T.right, h + (z ≤ T.Upper Limit? 1:0))

8:    else if z ≤ S then

9:       return Path Length(x, T.left, h + (z > T.Lower Limit? 1:0))

10:    end if

To assess the path length h effectively [39], the acceptable range of z in trained iTree is adopted to capture any unseen anomaly pixels that are out of range in the evaluation stage. This strategy utilizes multi-dimension spectral information (z) to make the DF more sensitive to anomaly pixels. As is shown in Figure 3, the upper and lower bounds are determined by the value of S and v in the training stage, where the value of v is the difference between the maximum and minimum of z in a sub-node. If a test pixel is outside the range in step 7 or step 10 of Algorithm 2, i.e., z ≤ Upper Limit or z > Lower Limit, the way of the path length without increasing punishes its anomaly score.

When a test pixel (x) reaches the exNode in step 1 of Algorithm 2, a total path length in each iTree is gained by the sum of h and c(T.Size), and the final average path length of t trees is formulated as

$$E\left(h\left(x\right)\right)=\frac{1}{t}{\displaystyle \sum _{i=1}^t{h}_i(x)},$$

(8)

where h_i(x) denotes the path length of a test pixel in the i-th tree of the Diforest. Finally, the anomaly score of E(h(x)) is normalized to (0, 1) by the formula of

$$s(x)=2^{-\frac{E\left(h(x)\right)}{c\left(n\right)}},$$

(9)

where c(n) is the normalizing factor that is calculated by the Equation (7).

3.3. Multiscale Fusion with Guided Filtering

As displayed in Figure 4a–c, although the detection results of the DF can capture the anomaly pixels well on the Gainesville dataset, a large number of background pixels are erroneously detected as the abnormal, resulting in the high FARs. To this end, the module of multiscale fusion with guided filtering is proposed in this paper, which can integrate effectively the multiscale DF detection maps with different spatial information to generate the final detection result. The implementation process of MFGF is shown in Figure 5.

First, the multiscale results from the DFs are normalized, and they are fused by pixel-level multiplication to obtain the result R₁. The fusion way can suppress the background, but a few anomaly pixels are missed. Then, the guided filtering [40] is applied to refine the details of the anomaly target, and it is defined as

$$O_i=a_k{I}_i+b_k,$$

(10)

where $O_i$ is a refined anomaly detection map, and the initial map of $O_i$ is R₁; $I_i$ is a guide map, i.e., the top three principal components $Y\in {ℝ}^{N\times 3}$; the $a_k$ and $b_k$ are linear coefficients, and they are defined as

$$\begin{array}{l}{a}_j={({\displaystyle \sum j}+\epsilon \mathrm{U})}^{-1}(\frac{1}{|\Omega |}{\displaystyle \sum _{i\in {\omega }_j}{I}_1^i}{\widehat{D}}_i-u_j{\tilde{D}}_j)\\ b_j={\tilde{D}}_j-a_j{u}_j\hspace{0.17em}\hspace{0.17em}{\widehat{D}}_i\end{array}$$

(11)

where $\epsilon$ means the regularization parameters, and it is set to 0.001; $u_j$ and ${\displaystyle \sum j}$ represent the mean matrix and covariance of guide image $I_i$ in the window ${\Omega }_j$ centered at the pixel j, respectively; $U$ is an identity matrix; $|\Omega |$ represents the number of pixels in the window ${\Omega }_j$; $\hat{D}$ denotes the refined anomaly detection map, i.e., R₁; ${\tilde{D}}_j$ is the mean of a refined map in the window ${\Omega }_j$, and the size of the window is set to $3\times 3$. More details are shown in [41,42]. Finally, the refined result O_i is combined with the initial fusion result of R₁ by pixel-level addition to get the result of R₂, and the final detection map R₃ is obtained by pixel-level multiplication of R₂ and O_i.

4. Experiment Results and Analysis

In this section, the experimental setup is introduced, and the qualitative and quantitative experiment comparisons are made to demonstrate the superiority of the proposed method. In addition, parameter analyses are carried out.

The experiments are implemented by a computer with Intel Core i7-9700 CPU at 3.00 GHz and 16 GB RAM (Intel, Santa Clara, CA, USA). The work platform of compared methods is MATLAB 2017b. The work platform of our method is python 3.6 with conda 4.5.4, and the applied packages are numpy 1.19.5, scipy 1.5.4, matplotlib 3.3.4, scikit-learn 0.24.2 and scikit-image 0.17.2.

4.1. Experiment Setup

4.1.1. Hyperspectral Dataset

To effectively verify the performance of the proposed method, we conduct experiments on four public hyperspectral datasets, i.e., Cat Island, Pavia, Gainesville, and San Diego. More details are described in

Table 1. Details of the Cat Island, Pavia, Gainesville, and San Diego datasets.

Dataset	Sensor	Spatial Resolution	Image Size	Anomaly Types
Cat Island	AVRIS	17.2 m/pixel	150 × 150 × 188	Airplanes
Pavia	ROSIS-03	1.3 m/pixel	150 × 150 × 102	Vehicles
Gainesville	AVRIS	3.5 m/pixel	100 × 100 × 191	Buildings
San Diego	AVRIS	3.5 m/pixel	120 × 120 × 189	Airplanes

.

Cat Island Dataset: the first dataset, which is captured by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) sensor on 12 September 2010. The dataset covers Gulfport, Southern Mississippi, USA. The size is 150 × 150 with 188 spectral bands, and the spatial resolution is 17.2 m/pixel. As shown in Figure 6(Ⅰ-a,Ⅰ-b), the anomaly object is an airplane with 19 pixels and accounts for 0.08% of the whole dataset.

Pavia Dataset: the second dataset is captured by the Reflective Optics System Imaging Spectrometer 03 (ROSIS-03) sensor, and its scene is the center of Pavia city in Northern Italy. The size is 150 × 150 with 102 spectral bands, and the spatial resolution is 1.3 m/pixel. The wavelength range is from 430 nm to 860 nm. Additionally, there are six anomaly targets in the dataset of Figure 6(Ⅱ-a,Ⅱ-b), and the number of anomaly pixels is 68 and 0.3% of the whole image.

Gainesville Dataset: the third dataset is captured by the AVIRIS sensor on 4 September 2010. The background is the city of Gainesville, FL, USA. The size of the dataset is 150 × 150, and the number of spectral bands is 102. The spatial resolution is 1.3 m/pixel. As is depicted in Figure 6(Ⅲ-a,Ⅲ-b), several buildings are anomaly objects, and it has 52 pixels and accounts for 0.52% of the whole dataset.

San Diego Dataset: the fourth dataset is captured by the AVRIS sensor, and the scene is the area of San Diego airport, CA, USA. The pixel size is 120 × 120, and the band number is 189. The spatial resolution is 3.5 m, and the range of wavelength is 400–2500 nm. Figure 6(Ⅳ-a,Ⅳ-b) are the pseudo color map and ground truth (GT), respectively. It is obvious that three airplanes are anomalies and have 57 pixels, which accounts for 0.4% of the total pixels of the dataset.

4.1.2. Compared Methods

The proposed algorithm is compared with six state-of-the-art algorithms in this paper, including RX [17], CRD [24], LRASR [9], LSDM-MoG [29], SSDF [15], and KIFD [13]. The RX is a classical statistics-based method. The CRD is a collaborative representation method. the LRASR and LSDM-MoG are low-rank representation methods. All settings of compared methods refer to the corresponding papers.

4.1.3. Evaluation Metrics

To effectively evaluate the performance of the proposed method, we apply four commonly used evaluation indicators, i.e., the area under the curve (AUC) of (P_d, P_f) and (P_f, τ), the receiver operating characteristic (ROC) curve of (P_d, P_f) and (P_f, τ), the time complexity and running time.

The ROC curve is to describe the tradeoff between the probability of detection (P_d) and the false alarm rate (P_f). When the threshold (τ) is set to an anomaly detection map, the probabilities of P_d and P_f can be obtained by the formula of

$$\begin{array}{l}{P}_d=\frac{N_D}{N_T}\\ P_f=\frac{N_F}{N}\end{array}$$

(12)

where $N_T$ denotes the total number of true target pixels in an HSI, and $N_D$ represents the number of true object pixels which are detected as objects; $N_F$ means the total number of background pixels in an HSI, and $N$ means denotes the total number of pixels in the HSI. By the effect of different threshold τ, the corresponding value of P_d and P_f can generate the ROC curve of (P_d, P_f) and (P_f, τ). Notably, the ROC curve of (P_d, P_f) is close to the upper left, and the performance is better; the ROC curve of (P_f, τ) is close to the lower left, and the detection effect is better.

The AUC, which is obtained by calculating the area of ROC, is applied to quantitatively evaluate the detection performance. Correspondingly, the value of AUC has (P_d, P_f) and (P_f, τ). When detection performance is better, the AUC of (P_d, P_f) is larger, and the AUC of (P_f, τ) is smaller.

In addition, the time complexity and running time are utilized to estimate the efficiency of different methods.

4.2. Detection Performance

In this section, the qualitative and quantitative comparisons are to evaluate the detection performance. For example, Figure 7 is the color detection map to make a quantitative comparison, and Table 2 and Table 3, Figure 8 and Figure 9 are the quantitative comparisons.

4.2.1. Qualitative Comparison

The qualitative comparison of detection color maps on four datasets is exhibited in Figure 7. Compared with other methods, the MSGDF method detects anomalous objects accurately and clearly in various complex backgrounds. For example, as for the Pavia detection maps of IF-based methods, the SSDF method cannot detect some anomaly objects of the buildings; in the KIFD method, although it can detect the buildings well, many background pixels are incorrectly detected; the proposed MSGDF method can detect all anomaly objects of buildings in the situation of low FARs. Hence, the proposed method has a better detection ability in subjective comparison.

Figure 7. Color detection maps on four datasets.

Figure 7. Color detection maps on four datasets.

4.2.2. Quantitative Comparison

The AUC comparisons are depicted in Table 2. From the AUC value in the single and entire dataset, the proposed MSGDF method has the best performance of the AUC of (P_d, P_f) and the AUC of (P_f, τ) in all methods. Meanwhile, compared with the algorithm of the second-best performance, the proposed algorithm has obvious advantages. Specifically, on the Cat Island dataset, the AUC score of (P_d, P_f) is 0.9992 and close to the ideal value of 1 in the proposed method, and our AUC score of (P_f, τ) is 0.0005 and close to 0, which is evidently superior to the CRD methods of the second-best performance. As for the challenging Pavia dataset, the RX has good detection performance, and its AUC value of (P_d, P_f) and (P_f, τ) are 0.9538 and 0.0233; in the proposed MSGDF method, our AUC value of (P_d, P_f) is 0.9868, which is 0.033 higher than the RX, and our AUC value of (P_f, τ) is 0.0014, which is 0.0216 lower than the RX. Furthermore, the AUC score of IF-based methods is listed in Table 2, we observe that the IF-based methods cannot acquire good performance of the AUC of (P_d, P_f) and (P_f, τ) at the same time. Taking the Gainesville dataset as an example, in the KIFD method, its AUC score of (P_d, P_f) and (P_f, τ) are 0.9921, and 0.1374, so we can think that the KIFD has a high AUC score of (P_d, P_f) and (P_d, P_f); on the contrary, the SSDF method has the low AUC score of (P_d, P_f) and (P_d, P_f), and its AUC score of (P_d, P_f) and (P_f, τ) are 0.8181 and 0.0506; in the proposed MSGDF method, it acquires a high AUC score of (P_d, P_f) and low AUC score of (P_d, P_f), the AUC score of (P_d, P_f) and (P_f, τ) are 0.9961 and 0.0047, respectively. In summary, the proposed method has a better detection rate and false alarm rate in the AUC comparison.

Table 2. AUC of (P_d, P_f) and (P_f,τ) on the four datasets. The bold mark and underlined mark indicate the first- and second-best performance, respectively.

Datasets	RX [17]	CRD [24]	LRASR [9]	LSDM-MoG [29]	SSDF [15]	KIFD [13]	MSGDF
AUC of (Pd, Pf)
Cat Island	0.9807	0.9916	0.4203	0.9942	0.9368	0.9897	0.9992
Pavia	0.9538	0.9537	0.9146	0.8646	0.8856	0.8164	0.9868
Gainesville	0.9513	0.9597	0.5897	0.9438	0.8181	0.9921	0.9961
San Diego	0.9111	0.9791	0.9853	0.9320	0.9884	0.9922	0.9960
Average value	0.9492	0.9710	0.7275	0.9337	0.9072	0.9476	0.9945
AUC of (Pf, τ)
Cat Island	0.0065	0.0058	0.3398	0.0181	0.0417	0.0952	0.0005
Pavia	0.0233	0.1104	0.1073	0.2708	0.0462	0.0815	0.0014
Gainesville	0.0351	0.0229	0.2979	0.0619	0.0506	0.1374	0.0047
San Diego	0.0406	0.0360	0.0975	0.0882	0.0399	0.0980	0.0021
Average value	0.0264	0.0438	0.2106	0.1098	0.0446	0.1030	0.0022

Table 2. AUC of (P_d, P_f) and (P_f,τ) on the four datasets. The bold mark and underlined mark indicate the first- and second-best performance, respectively.

Datasets	RX [17]	CRD [24]	LRASR [9]	LSDM-MoG [29]	SSDF [15]	KIFD [13]	MSGDF
AUC of (Pd, Pf)
Cat Island	0.9807	0.9916	0.4203	0.9942	0.9368	0.9897	0.9992
Pavia	0.9538	0.9537	0.9146	0.8646	0.8856	0.8164	0.9868
Gainesville	0.9513	0.9597	0.5897	0.9438	0.8181	0.9921	0.9961
San Diego	0.9111	0.9791	0.9853	0.9320	0.9884	0.9922	0.9960
Average value	0.9492	0.9710	0.7275	0.9337	0.9072	0.9476	0.9945
AUC of (Pf, τ)
Cat Island	0.0065	0.0058	0.3398	0.0181	0.0417	0.0952	0.0005
Pavia	0.0233	0.1104	0.1073	0.2708	0.0462	0.0815	0.0014
Gainesville	0.0351	0.0229	0.2979	0.0619	0.0506	0.1374	0.0047
San Diego	0.0406	0.0360	0.0975	0.0882	0.0399	0.0980	0.0021
Average value	0.0264	0.0438	0.2106	0.1098	0.0446	0.1030	0.0022

The ROC curves of (P_d, P_f) and (P_f, τ) are described in Figure 8 and Figure 9. As for the ROC curves of (P_d, P_f), we can observe that the curve of the proposed MSGDF algorithm is always at the top left of six compared algorithms on Cat Island, Pavia, and Gainesville datasets. In the San Diego dataset, A small part on the left of our curve is below the curve of the KIFD and LRASR methods, and the rest part of our curve is above the compared methods. From the overall curve, our method is also at the top left of the other methods. Therefore, the ROC curves of (P_d, P_f) illustrate that the proposed MSGDF method is superior to the compared methods. In the dimension of the ROC curves of (P_f, τ), it is evident that the curve of the proposed MSGDF is at the bottom left of six state-of-the-art algorithms, which indicates that our method has a better performance on the false alarm rates.

These experiment results illustrate that the proposed method has preeminent detection accuracy and FARs, which demonstrates that the design of MSGDF can better dig out HSI spectral–spatial information to detect local anomalies and improve the problem of high FARs.

Figure 8. ROC curves of (P_d, P_f) on the four datasets. (a) Cat Island, (b) Pavia, (c) Gainesville, (d) San Diego.

Figure 8. ROC curves of (P_d, P_f) on the four datasets. (a) Cat Island, (b) Pavia, (c) Gainesville, (d) San Diego.

Figure 9. ROC curves of (P_f, τ) on the four datasets. (a) Cat Island, (b) Pavia, (c) Gainesville, (d) San Diego.

Figure 9. ROC curves of (P_f, τ) on the four datasets. (a) Cat Island, (b) Pavia, (c) Gainesville, (d) San Diego.

The time complexity analysis of the proposed DF model is carried out, and the proposed DF model mainly consists of three parts, i.e., calculating the hyperplane projection value, computing the gain splitting criterion, and partitioning the projection attribute value. Therefore, the time complexity is O(tn(qn + logn + n)) in the training stage, and the time complexity of the evaluation stage is O(tnqN). The total time complexity of the proposed DF model is O(tn(q(n + N) + logn + n)). In addition, the comparison results of the running time on four HSI datasets are listed in Table 3, and it can be observed that the proposed MSGDF method is superior to the other methods, except the RX method. In the comparison of the IF-based methods, our method is the fastest. This is because the overall design of the proposed MSGDF model is relatively simple and valid to reduce time consumption. Two other IF-based methods are relatively complex, causing increasing the running time. Concretely, the KIFD algorithm adopts the kernel space and the iteration strategy to optimize its detector; the SSDF algorithm applies the axis-parallel subspace selection process and complex splitting criterion to enhance the discriminant ability. The strategies of KIFD and SSDF methods trade the lower running efficiency for high detection accuracy.

Table 3. Running time on the four datasets. The bold mark and underlined mark indicate the first- and second-best performance, respectively.

Datasets	RX	CRD	LRASR	LSDM-MoG	SSDF	KIFD	MSGDF
Cat Island	0.9065	16.4834	2.9095	47.4954	47.4678	71.2541	15.3171
Pavia	0.1729	14.7069	92.2870	20.2075	58.2189	87.1228	14.0285
Gainesville	0.4846	3.1657	1.5270	19.9304	40.4361	64.1250	3.3587
San Diego	0.5166	15.0573	1.8886	13.0941	47.2768	68.3000	6.3185
Average value	0.52015	12.35333	24.65303	25.18185	48.3499	72.70048	9.7557

Table 3. Running time on the four datasets. The bold mark and underlined mark indicate the first- and second-best performance, respectively.

Datasets	RX	CRD	LRASR	LSDM-MoG	SSDF	KIFD	MSGDF
Cat Island	0.9065	16.4834	2.9095	47.4954	47.4678	71.2541	15.3171
Pavia	0.1729	14.7069	92.2870	20.2075	58.2189	87.1228	14.0285
Gainesville	0.4846	3.1657	1.5270	19.9304	40.4361	64.1250	3.3587
San Diego	0.5166	15.0573	1.8886	13.0941	47.2768	68.3000	6.3185
Average value	0.52015	12.35333	24.65303	25.18185	48.3499	72.70048	9.7557

4.3. Parameter Analysis

4.3.1. Parameter Analysis of Multiscale Superpixel Segmentation

The parameter setting in the number of superpixels is crucial for the HAD because it determines the feature extraction of HSI spatial information. In the article, the number of superpixels is set from 1 to 19, and the reason is that the range not only effectively distinguishes different types of homogeneous scene regions but also ensures the subsample number of each DF [14]. For the number of superpixel scale, its range is from 2 to 6, because this range can significantly reduce the false alarm rate in this paper. Moreover, to obtain the setting of multiscale superpixel segmentation, the single scale superpixel setting is acquired by selecting the best detection performance, and then based on the single scale superpixel setting, the number of multiscale superpixels is obtained by the best detection performance.

The AUC score of (P_d, P_f) of the single-scale superpixel DF is displayed in Figure 10. It can be observed that the superpixel-guided DF in the situation that the number of superpixels is over 1 is better than the original DF in the situation that the number of superpixels is 1 in the overall detection effect. Especially, when the optimum number of superpixels is selected, the advantage of single-scale superpixel DF is more obvious. This experiment results verify the necessity of the superpixel-guided DF in HAD. According to the best detection accuracy of (P_d, P_f) on four HSI datasets, the number of superpixels on Cat Island, Pavia, Gainesville, Los Angeles, and San Diego datasets are set to 17, 3, 9, and 5, respectively. Furthermore, we find that the best number setting of single-scale superpixels on each HSI data set is affected by the HSI scene complexity.

To facilitate the setting of multiscale superpixels on four HSI datasets, the superpixel number of the obtained single-scale setting is the center, which can set the other numbers of the multiscale superpixels. For the superpixel setting of every scale, we select the best detection accuracy in the superpixel number combination as the scale setting result. On the basis of the principles, the AUC scores of the multiscale superpixel setting (MSS) are listed in

Table 4. Number of superpixel scales for the influence of AUC in MSGDF. The bold mark indicates the best performance.

Scale		2	3	4	5	6
	MSS	(15, 17)	(15, 17, 19)	(13–19)	(11–19)	(9–19)
Cat Island	AUC of (Pd, Pf)	0.9976	0.9992	0.9990	0.9990	0.9992
	AUC of (Pf, τ)	0.0033	0.0005	0.0001	0.0001	0
	MSS	(3, 5)	(1, 3, 5)	(7–13)	(7–15)	(7–17)
Pavia	AUC of (Pd, Pf)	0.9807	0.9961	0.9958	0.9959	0.996
	AUC of (Pf, τ)	0.0065	0.0047	0.0041	0.0016	0.001
	MSS	(9, 11)	(9, 11, 13)	(7–13)	(7–15)	(7–17)
Gainesville	AUC of (Pd, Pf)	0.9940	0.9961	0.9958	0.9959	0.996
	AUC of (Pf, τ)	0.0224	0.0047	0.0041	0.0016	0.001
	MSS	(5, 7)	(3, 5, 7)	(1–7)	(1–9)	(1–11)
San Diego	AUC of (Pd, Pf)	0.9948	0.996	0.9948	0.9946	0.9942
	AUC of (Pf, τ)	0.0087	0.0021	0.0024	0.0008	0.0004
Average Value	AUC of (Pd, Pf)	0.9918	0.9945	0.9930	0.9929	0.9930
	AUC of (Pf, τ)	0.0102	0.0022	0.0018	0.0007	0.0004

. From the AUC value of (P_f, τ), we can find that the false alarm rate decreases with the increase of the superpixel scales, and the best performance of FARs is got in the scale setting of 6. In the AUC score of (P_d, P_f), the detection accuracy of superpixel scale 3 is higher than the other scale. Finally, by considering comprehensively the AUC scores of (P_d, P_f) and (P_f, τ), when the number of the superpixel scale is 3, it can acquire the best performance of detection rates in the situation of low FARs. Therefore, the number settings of multiscale superpixels on the Cat Island, Pavia, Gainesville, Los Angeles, and San Diego datasets are set to (15, 17, 19), (1, 3, 5), (9, 11, 13), and (3, 5, 7), respectively.

4.3.2. Parameter Analysis of Discrimination Forest

The effect of three parameter settings in DF is analyzed in Figure 11, including the dimension (q), subsampling size (n), and the number of trees (t). To display clearly the impact of parameters, the red curve of average value is to represent the holistic detection impact, and the curve of other colors is to depict the changing trend in a single dataset. Significantly, the larger the value of the three parameters are, the more running time will be consumed.

Figure 11a depicts the influence of the subsample size n on AUC of (P_d, P_f), and the subsample setting range of DF is from 3% to 50% in each superpixel block. The AUC (P_d, P_f) value on the Cat Island dataset remains around 0.99; regarding the other three datasets, the AUC value begins to rise, then decreases, and finally tends to remain unchanged. Most notably, when the subsample is about 10–20%, the best performance of AUC (P_d, P_f) is gained on four datasets.

The number of dimensions q for the hyperplane building of DF is from 1 to 100. As is depicted in Figure 11b, the AUC values of (P_d, P_f) initially increase and then fluctuate in a small range on the four datasets, and the max AUC value of (P_d, P_f) is achieved at dimension number of 2. Furthermore, the changing trend of AUC curves illustrates that utilizing multi-dimension spectral bands is better than using a single-dimension spectral band in the DF. This phenomenon demonstrates the effectiveness of multi-dimension design on the gain split criterion.

The value of the number of discriminative trees is from 10 to 500 in Figure 11c. We observe that the changing trend of the AUC value of (P_d, P_f) is similar to Figure 11b, and the difference is that the AUC value of the Cat Island dataset changes slightly at 0.999.

Finally, according to comprehensive consideration of the average detection accuracy and time consumption, the parameters of dimension and trees are set to 2 and 100 respectively, and the subsample size is 10% of pixels in each HSI superpixel block.

5. Discussion

To illustrate the effectiveness of DF, the IF and proposed DF are tested on four HSI datasets, and their results with the AUC of (P_d, P_f) are displayed in Figure 12. From the detection accuracy results, the proposed DF is superior to the original IF both in the AUC (P_d, P_f) of the single dataset and average AUC (P_d, P_f) of the entire dataset, and the average AUC value of the DF method is 2% higher than the average AUC value of IF method. This verifies that the design of the DF is more sensitive to anomaly pixels than the design of the original IF. Meanwhile, it illustrates that the strategy of the gain split criterion and acceptable range is effective in the DF model.

To further explain the influence of improved strategies in our proposed DF, such as the gain split criterion and acceptable range, the experiment of component analysis is conducted in four hyperspectral data sets, and the results of detection accuracy are listed in

Table 5. The proposed DF component analysis by AUC of (P_d, P_f).

Optimization Strategies		Hyperspectral Data Set				Average
Gain Split Criterion	Acceptable Range	Cat Island	Pavia	Gainesville	San Diego
		0.9685	0.9409	0.9596	0.9747	0.9609
√		0.9800	0.9668	0.9680	0.9798	0.9737
	√	0.9837	0.9608	0.9610	0.9821	0.9719
√	√	0.9855	0.9760	0.9773	0.9832	0.9805

. Under the strategy of gain split criterion or acceptable range, the average AUC score of (P_d, P_f) of improved DF is over 1% higher than without improved DF, meanwhile, the AUC score rises up on each hyperspectral data set. These results illustrate the efficiency of two design strategies. On the other side, compared with adopting a single optimization strategy, the way of gain split criterion combined with acceptable range further improves AUC accuracy in each data set and the entire data sets, which verifies the rationality of combination of two optimization strategies.

6. Conclusions

A novel detection strategy with the MSGDF model is proposed for HAD. In this model, the multiscale superpixel segmentation is applied to excavate spatial features and guide local anomaly detection, and the split criterion of multi-dimension spectral bands is adopted to optimize each node division for the DF. The strategy, which combines DF with SLIC, makes full use of HSI spectral and spatial information to solve the problem of local anomalies. In addition, the MFGF is put forward to effectively reduce FARs of initial DF results. Finally, the experimental analysis validates the rationality of the MSGDF method, and the detection results illustrate that the MSGDF has better detection ability than state-of-the-art algorithms. Although the performance of our proposed algorithm is satisfactory, some meaningful problems remain to be explored in our future work. Specifically, the obtained background pixels are constructed semi-supervised iforest to enhance the sensitivity of anomalies in test stage. Furthermore, whether the IF methods can combine with the other HAD methods to improve detection effect, such as LRR, CR and AE.