Spatial–Spectral Joint Hyperspectral Anomaly Detection Based on a Two-Branch 3D Convolutional Autoencoder and Spatial Filtering

Abstract

Hyperspectral anomaly detection (HAD) is an important application of hyperspectral images (HSI) that can distinguish anomalies from background in an unsupervised manner. As a common unsupervised network in deep learning, autoencoders (AE) have been widely used in HAD and can highlight anomalies by reconstructing the background. This study proposed a novel spatial–spectral joint HAD method based on a two-branch 3D convolutional autoencoder and spatial filtering. We used the two-branch 3D convolutional autoencoder to fully extract the spatial–spectral joint features and spectral interband features of HSI. In addition, we used a morphological filter and a total variance curvature filter for spatial detection. Currently, most of the datasets used to validate the performance of HAD methods are airborne HSI, and there are few available satellite-borne HSI. For this reason, we constructed a dataset of satellite-borne HSI based on the GF-5 satellite for experimental validation of our anomaly detection method. The experimental results for the airborne and satellite-borne HSI demonstrated the superior performance of the proposed method compared with six state-of-the-art methods. The area under the curve (AUC) values of our proposed method on different HSI reached above 0.9, which is higher than those of the other methods.

1. Introduction

Hyperspectral imaging technology has the unique advantages of a high spectral resolution and a high spatial resolution, and it is one of the most significant scientific and technological breakthroughs since the development of remote sensing technology [1]. Hyperspectral images (HSI) are image cubes, enabling the acquisition of fine spectral information along with spatial information about the feature [2,3]. HSI has been widely used in different fields of remote sensing, such as classification [4], target detection [5], mineral resources surveys [6], gas emissions surveys [7], and anomaly detection [8]. Anomalies in HSI can be considered as pixels that are different from the background and usually have the following characteristics: (1) they have significant spectral differences from the surrounding background; (2) they are spatially isolated and often occupy fewer pixels or sub-pixels; (3) the spatial brightness of anomalous pixels in different bands differs from the surrounding background, being either brighter or darker [9,10,11]. Compared with target detection, which requires a priori information, anomaly detection can detect an anomaly without using any a priori information from the HSI, providing important information for subsequent image analysis. Therefore, hyperspectral anomaly detection, which is more practical [12,13], has been widely applied in many fields of remote sensing, such as target detection, mineral exploration and environmental monitoring [1,14]. In the past two decades, researchers have proposed many effective hyperspectral anomaly detection methods from different perspectives.

The statistical theory-based hyperspectral anomaly detection methods are the most classical, which use the statistical features of hyperspectral images to detect anomalies. The best-known and most representative statistical-based method is the Reed–Xiaoli (RX) algorithm [15]. The RX algorithm assumes that the background of the hyperspectral image follows a Gaussian distribution and obtains the anomaly detection result by calculating the Mahalanobis distance between the mean spectrum and the test pixel’s spectrum. Many improved methods based on the RX algorithm have been proposed. Multiple-window RXD (MW-RXD) was proposed to obtain better detection results when the anomalies have different shapes and sizes [16]. To address the problem that real hyperspectral images do not always follow a Gaussian distribution, Matteoli proposed a data-driven strategy for automatically estimating the background probability density function [17]. To reduce the effect of anomalous pixels on the estimation of background statistical features, weighted-RXD (W-RXD) sets different weights for background pixels and anomaly pixels; linear filter-based RXD (LF-RXD) filters out anomalies by calculating the possibility of each pixel being a background sample to [18]. RSAD [19] randomly selects background pixels several times to calculate the background statistics for anomaly detection. In addition, the kernel RX (KRX) algorithm [20] maps the raw data nonlinearly to a high-dimensional feature space where the background and anomalies can be separated more effectively. For real hyperspectral images, accurately modelling the background is statistically difficult, so many approaches from other perspectives have been proposed.

The distance-based detection method groups all pixels in a hyperspectral image according to their distance, and pixels that deviate from the centre of the group are considered to be anomalous pixels. The support vector description (SVDD) anomaly detection method calculates the amount of support vector description for each pixel and considers pixels outside the support region as anomalous pixels [21]. An anomaly detection method based on the local joint subspace and the support vector machine (LJSSVM) [22] was proposed by combining the SVM into a statistical approach. Clustering-based anomaly detection (CBAD) [23] uses a clustering algorithm to divide the original hyperspectral image into clusters, with a pixel further away from the centre of the nearest cluster considered to be more anomalous. The authors of [24] applied graph theory to the detection of hyperspectral anomalies.

The representation-based hyperspectral anomaly detection methods have become a hot research topic. The basic idea is that background pixels can be represented with a small level of error using dictionary atoms, while anomalous pixels cannot. The local sparsity divergence anomaly detection method (LSDAD) [25] provides a consistent sparse divergence index (SDI) and fuses local spectral sparse divergence with local spatial sparse divergence to represent the anomaly degree of the pixel. Ling et al. added sum-to-one and non-negative constraints to the abundance vector on top of the sparse representation model to ensure that it was physically meaningful [26]. Unlike sparse representation methods, the collaborative representation anomaly detection (CRD) method [27] represents the central pixel collaboratively with all pixels in the spatial neighbourhood and enhances the collaboration of interpixel representation with an L2 norm constraint so that every pixel participates in the representation.

Many matrix decomposition-based anomaly detection methods have been proposed due to the global low rank of the background of the hyperspectral image and the global sparsity of the anomaly. A low-rank and sparse matrix decomposition detection method (LRaSMD) [28] was proposed, which uses the Go Decomposition (GoDec) algorithm to solve the low-rank matrix decomposition problem and obtains the result by calculating the Euclidean distance of sparse matrix. Based on the LRaSMD method, the LSMAD method was proposed, which calculates the Mahalanobis distance using the background matrix obtained from the decomposition to obtain detection results [29]. Considering that the backgrounds of hyperspectral images belong to different subspaces, many methods based on the low-rank representation (LRR) have been proposed. To preserve the local geometric structure and the spatial relationships of HSI, an anomaly detection method based on graphs and total variational regularisation LRR (GTVLRR) has been proposed [30]. To give a more accurate representation of the background, the low-rank and sparse representation method (LRASR) adds a sparse regularisation term to the representation of the coefficient matrix based on the low-rank constraint on the background matrix and uses the K-means algorithm to construct the dictionary [31].

Deep learning-based methods are being increasingly applied to hyperspectral anomaly detection due to they can capture high-dimensional, nonlinear features and fit complex functions [32]. Du et al. applied a convolutional neural network (CNN) to detect hyperspectral anomalies and proposed an anomaly detection framework based on the migrating depth of convolutional neural networks to measure the degree of anomaly as the similarity between the pixel pair consisting of the pixel to be measured and the surrounding pixels [33]. Song et al. used CNN to obtain abundance maps to use as input for LRR to complete anomaly detection [34]. However, the lack of labeled anomalous samples in anomaly detection tasks imposes significant limitations on supervised methods. As one of the unsupervised neural networks, the autoencoder (AE) is increasingly being used for the detection of hyperspectral anomalies due to its not requiring labeled anomalous samples. To preserve the geometric structure and the local spatial consistency of HSI simultaneously, a method named robust graph autoencoder (RGAE) has been proposed [35]. The spectral constraint adversarial autoencoder approach (SC-AAE) [36] incorporates a spectral constraint strategy into adversarial autoencoders to fully utilize the spectral information of the hyperspectral data to extract the features of high-dimensional spectral vectors. To attenuate the effects of noise, interband nonlinear correlations, and other factors on detection, a stacked denoising autoencoder-based detection method (HADSDA) was proposed [37]. A number of autoencoder-based methods of detecting hyperspectral anomalies are now available and have had good results. However, they still have some problems. Firstly, flattening the hyperspectral images’ 3D cube data into 2D matrix data results in the loss of information on the spatial–spectral structure. Secondly, the rich spectral interband information of HSI is underutilized. Finally, the presence of anomalous pixels in the training set will make the reconstruction error for anomalies smaller, resulting in less ability to distinguish between the background and an anomaly.

As the AE network can learn the background feature well and reconstruct the background without using labeled anomalous samples, it is used as the basic model in this paper for hyperspectral anomaly detection. In order to address the abovementioned issues, we proposed a spatial–spectral joint HAD method based on a two-branch 3D convolutional autoencoder(3D-CAE) and spatial filtering. In addition, to address the lack of satellite-borne HSI for anomaly detection experiments, we constructed a satellite-borne hyperspectral dataset containing HSI with different backgrounds and different anomalous targets. The main innovative contributions of our work can be summarised as follows.

(1)

A novel two-branch 3D-CAE was developed to fully extract the spatial–spectral joint features and spectral interband features of HSI, and novel multi-scale spectral difference data were used as the input of the second network branch.

(2)

A morphological filter and a total variance curvature filter were used for spatial detection, and the spatial detection result was also used to filter the background sample set for training the network.

(3)

A satellite-borne hyperspectral dataset based on the images acquired by the GF-5 satellite was constructed that can be used to validate the effectiveness of many HAD methods. We used six state-of-the-art methods to demonstrate the validity of the proposed method, not only with the commonly used airborne hyperspectral images but also with satellite-borne HSI.

The rest of the study is structured as follows. Section 2 details the two parts of the hyperspectral anomaly detection method proposed in this study. In Section 3, we validate the effectiveness of the proposed method using real airborne and satellite-borne HSI, and Section 4 presents the conclusions.

2. Proposed Method

This section describes in detail the method of detecting hyperspectral anomalies proposed in this study. As shown in Figure 1, the proposed method consists of two parts: the spatial anomaly detection part and the spatial–spectral joint detection part. First, in the spatial anomaly detection part, morphological filtering and total variance curvature filtering [38] are applied to the first few principal components of the original HSI to extract the spatial features of anomalies and suppress the background information, thus obtaining the spatial detection results. Then, in the spatial–spectral joint detection part, we constructed a two-branch 3D convolutional autoencoder network [39]. To better extract the spectral and interspectral information of the hyperspectral image, we proposed a multi-scale spectral difference feature as the input for the second branch of the network. The original HSI serves as the input of the first branch. After the training was completed, we input the original HSI and the multi-scale spectral difference data into the network and took the reconstruction error of the first branch as the results of spatial–spectral joint anomaly detection result. Finally, we fused the results of the two parts to obtain the final result of anomaly detection.

2.1. Spatial Detection

Based on the characteristics of anomalies described earlier, we know that anomalies typically have rich spatial information that can be utilized. To make spatial information more significant, we first performed principal component analysis (PCA) on the original hyperspectral images to obtain the first few principal components. The PCA method can concentrate the main spatial information in the first few principal components [40]. We then applied a morphological gradient and top-hat operations [41] to the first few principal components to obtain the boundaries of the anomalies and isolated pixels (i.e., potentially anomaly pixels) with higher brightness values than the surrounding area in the first few components.

In this study, the original hyperspectral image is represented as $X\in {ℝ}^{\mathrm{M}\times \mathrm{N}\times \mathrm{L}}$, where $\mathrm{M}$ and $\mathrm{N}$ are the sizes of the spatial dimensions of the HSI, and $\mathrm{L}$ is the number of bands in the spectral dimension. The first few principal components $C_i\in {ℝ}^{\mathrm{M}\times \mathrm{N}\times 1}$ were obtained by using PCA for $X$:

$$C_i={\mathrm{PC}}_i\left(X\right),i=1,2\cdots \mathrm{P}$$

(1)

where ${\mathrm{PC}}_i\left(X\right)\in \left\{{\mathrm{PC}}_1\left(X\right),{\mathrm{PC}}_2\left(X\right),\cdots ,{\mathrm{PC}}_{\mathrm{B}}\left(X\right)\right\}$ are the principal components of the hyperspectral image obtained using PCA. Then, morphological gradient and top-hat operations were applied to $C_i$:

$$G_i=grad\left(C_{\mathit{i}}\right),T_{\mathit{i}}=tophat\left(C_{\mathbf{i}}\right)$$

(2)

$$grad\left(C_i\right)=\left(C_{\mathit{i}}\oplus \mathrm{SE}\right)-\left(C_{\mathit{i}}\ominus \mathrm{SE}\right)$$

(3)

$$tophat\left(C_i\right)=C_{\mathit{i}}-\left(C_{\mathit{i}}\ominus \mathrm{SE}\right)\oplus \mathrm{SE}$$

(4)

where $grad\left(·\right)$ and $tophat\left(·\right)$ are the morphological gradient and top-hat operations, $\oplus$ denotes the morphological expansion operation, $\ominus$ denotes the morphological erosion operation, and $\mathrm{SE}$ denotes the structural element (we used the cross-structure element). After the morphological feature map had been obtained, the morphological feature map was normalized and fused to obtain the initial result of spatial detection $I\in {ℝ}^{\mathrm{M}\times \mathrm{N}}$. The value of $I$ can be calculated as

$$I=\frac{1}{\mathrm{P}}{\displaystyle \sum }_{i=1}^{\mathrm{P}}\frac{\left(G_{\mathrm{i}}+T_{\mathrm{i}}\right)}{\mathbf{2}}$$

(5)

The operation above obtained the original result of spatial anomaly detection while also obtaining some background information. Therefore, the background information needed to be further subtracted from $I$ to obtain the final result of spatial detection $A_1$. The total variance curvature filter (TVCF) [38,42] assumed that the surface of an image is a piecewise constancy surface and that the HSI has inherent piecewise constancy due to their sparsity. Therefore, in this study, the background image, denoted as $B$, was obtained by multiple iterations of filtering using the TVCF for $I$.

In the TVCF, a domain decomposition method is used to divide all pixels of the image into four subsets such that neighboring pixels in a four-connected neighborhood belong to different subsets. The specific decomposition method is as follows: first, all image pixels are divided into two non-adjacent subsets: the “white” points W and the “black” points B. Then, each of these two subsets is further divided into triangles and circles. This results in four subsets: black circles B_C, black triangles B_T, white circles W_C and white triangles W_T, as shown in Figure 2. Then, for all pixels in each subset, the projection of the pixel to the eight surfaces within its 3 × 3 neighborhood was calculated. The projection with the smallest absolute value was then selected from these projections and added to the original pixel to obtain the updated result. The algorithm is summarized in detail in Algorithm 1. After several iterations using TVCF, the background image $B$ could be obtained, subtracted from the initial spatial anomaly feature map $I$, and the difference was squared to enhance the attenuating background of the anomalies to obtain the results of spatial detection.

$$A_1={\left(I-B\right)}^2$$

(6)

Algorithm 1. Total variation curvature filter

Input: $I_{m,n}\in B_C,B_T,W_C,W_T$

1: ${\widehat{d}}_1=(I_{m-1,n-1}+I_{m-1,n}+I_{m,n-1}+I_{m+1,n-1}+I_{m+1,n})/5-I_{m,n}$

2: ${\widehat{d}}_2=(I_{m-1,n}+I_{m-1,n+1}+I_{m,n+1}+I_{m+1,n}+I_{m+1,n+1})/5-I_{m,n}$

3: ${\widehat{d}}_3=(I_{m-1,n-1}+I_{m-1,n}+I_{m-1,n+1}+I_{m,n-1}+I_{m,n+1})/5-I_{m,n}$

4: ${\widehat{d}}_4=(I_{m+1,n-1}+I_{m+1,n}+I_{m+1,n+1}+I_{m,n-1}+I_{m,n+1})/5-I_{m,n}$

5: ${\widehat{d}}_5=(I_{m-1,n-1}+I_{m-1,n}+I_{m-1,n+1}+I_{m,n-1}+I_{m,n+1})/5-I_{m,n}$

6: ${\widehat{d}}_6=(I_{m-1,n-1}+I_{m-1,n}+I_{m-1,n+1}+I_{m,n-1}+I_{m+1,n+1})/5-I_{m,n}$

7: ${\widehat{d}}_7=(I_{m+1,n-1}+I_{m+1,n}+I_{m+1,n+1}+I_{m-1,n-1}+I_{m,n-1})/5-I_{m,n}$

8: ${\widehat{d}}_8=(I_{m+1,n-1}+I_{m+1,n}+I_{m+1,n+1}+I_{m-1,n+1}+I_{m,n+1})/5-I_{m,n}$

9: Find $d_m$$, \mathrm{such} \mathrm{that} \left|d_m\right|=\min \left\{\left|{\widehat{d}}_i\right|,i=1,\cdots ,8\right\}$

Output: $I_{m,n}=I_{m,n}+d_m$

2.2. Spatial–Spectral Joint Detection

AE networks can use spectral reconstruction errors to measure the degree of anomaly of a pixel. However, the traditional AE networks flatten the 3D cube data into a 2D matrix and use 1D spectral vectors as inputs to the network, ignoring the spatial and joint spatial–spectral information. For spatial–spectral joint detection, we proposed a two-branch 3D-CAE network and a multi-scale spectral difference feature in order to make full use of the spatial–spectral joint information. To avoid the involvement of anomalous pixels in the training of the network, we proposed a coarse screening strategy for the background patch samples. Below, we present a detailed description of the multi-scale spectral differencing and network inputs, the network’s structure, the loss function, and the final results of anomaly detection.

2.2.1. Multi-Scale Spectral Difference Feature Data and Network Inputs

Spectral differencing is a commonly used strategy for analyzing HSI’s spectral features, where the spectral differences at different scales represent different spectral features. Smaller spectral differences can reflect the details of the spectral changes, while larger spectral differences can reflect the trends of the spectral changes. The difference between an anomaly and the background may be greater at some scales. Figure 3 shows the different scales of the spectral differences in the data. In order to fully extract the spectral features at different scales, we proposed a multi-scale spectral difference strategy, where spectral difference data at different scales are calculated for the original hyperspectral image and fused to obtain multi-scale spectral difference data as the input for one branch of the two-branch 3D-CAE network.

The features of spectral difference at different scales can be expressed as

$$s_i\left(:,:\right)={\displaystyle \sum }_{k=1}^{L-i}\left|X\left(:,:,k+i\right)-X\left(:,:,k\right)\right|$$

(7)

where $i\in \left[1,L-1\right]$ denotes the scale and $s_i\left(:,:\right)\in {ℝ}^{\mathrm{M}\times \mathrm{N}}$ is the spectral difference at scale $i$. Then, the multi-scale spectral difference data cube can be obtained by combining the spectral difference at different scales, described as

$$S\left(:,:,i\right)=s_i\left(:,:\right)$$

(8)

where $S\in {ℝ}^{M\times N\times L-1}$ and, originally, the HSI data $X\in {ℝ}^{\mathrm{M}\times \mathrm{N}\times \mathrm{L}}$ are the two input datasets for the two branches of network.

Further, to avoid patches containing abnormal pixels being fed into the network, patches in $X$ and $S$ were not fed directly into the network. We used the spatial detection result $A_1$ to screen the background patches in $X$ and $S$ within the network. First, the background anomaly segmentation map $A_{1\_seg}$ was obtained using the Otsu thresholding method [43] for $A_1$. Then, the patches in $X$ and $S$ whose center pixels corresponded to a value of zero in $A_{1\_seg}$ were considered to be background patches and were inputted into the network. Thus, the inputs of the two branches of the network were $x_{m,n}$ and $s_{m,n}$, which are described as

$$x_{m,n}\in X=\left\{x_{m,n}|A_{1\_seg}\left(m,n\right)=0\right\}$$

(9)

$$s_{m,n}\in S=\left\{s_{m,n}|A_{1\_seg}\left(m,n\right)=0\right\}$$

(10)

where $x_{m,n},s_{m,n}$ are the cube patch in $X$ and $S$ centered at $\left(m,n\right)$ with sizes of $p\times p\times L$ and $p\times p\times L-1$, respectively.

2.2.2. Architecture of the Two-Branch 3D-CAE Network

Hyperspectral images are three-dimensional cube data, while the AE network flattens the cube data into a two-dimensional matrix during training, which can only extract spectral features and lose the spatial-spectral joint features of hyperspectral images. The 3D-CAE architecture is a variant of the traditional AE network, which takes the three-dimensional cube data as input and can extract both spectral and spatial features simultaneously. Therefore, 3D-CAE is suitable for hyperspectral-related applications and has been successfully applied in hyperspectral image processing, such as classification and detection [39,44,45,46]. To exploit both the rich spectral features and the spatial–spectral joint features of HSI, we were inspired by the research on 3D-CAE and proposed a two-branch 3D-CAE network. The network’s architecture is shown in Figure 4. The first branch of the network extracts the deep spatial–spectral joint features from the original HSI $X$, while the second branch extracts the spectral features from the multi-scale spectral difference data $S$. The spectral difference branch improves the ability of the model to distinguish between different spectra. For the two different branches, the network’s architecture is the same. The specific configuration of the model is as follows.

In the encoding part, the features of two patches $x_{m,n}$ and $s_{m,n}$ were extracted separately from the two branches. The input 3D patches of the two branches of the network were $x_{m,n}$ and $s_{m,n}$, with a size of $p\times p\times L$ and $p\times p\times L-1$. We set the size a $p=5$ for both $x_{m,n}$ and $s_{m,n}$ because anomalies typically have a smaller spatial size (5 × 5 window is sufficient to encompass the anomalies), and it is crucial to ensure that two branches extract features at the same spatial scale. We used two 3D convolutional layers to extract the spatial–spectral joint features. The first convolutional layer had 32 convolutional kernels with a size of $3\times 3\times 5$ and a stride of $1\times 1\times 1$. The convolutional layer was followed by the batch normalization layer and the PRelu layer. The second convolutional layer had 64 convolutional kernels with a size of $3\times 3\times 5$ and a stride of $1\times 1\times 1$. This convolutional layer was followed by the batch normalization layer and a sigmoid layer. The features obtained from the two convolutional layers were flattened to a 1D vector which became the input of the fully connected layer. The output of this layer of the two branches was concatenated and passed through a fully connected layer again to obtain the latent vector.

In the decoding part, the latent layer vector was first restored to its original scale by a fully connected layer. It was then separated and transformed back into two 3D patches, which were used as input for the decoding part of the two branches. The architecture of the network of the decoding part was a mirror image of the encoding part, except that a transposed convolutional layer was used instead of a convolutional layer. The output of the network ${\widehat{x}}_{m,n}$ and ${\widehat{s}}_{m,n}$, which had the same size as $x_{m,n}$ and $s_{m,n}$, was the reconstructed 3D patch.

2.2.3. Loss Function and Final Results of Anomaly Detection

The loss functions of the two branches were calculated the same way and were summed to obtain the loss function of the network. With the original image branch as an example, the loss function had two components: $L_{spatial}$, which measures the input patch’s spatial similarity to the reconstructed patch, and $L_{spectral}$, which measures the input patch’s spectral similarity to the reconstructed patch. Since the background screening session made the central pixel of each patch more likely to be the background, we multiplied $L_{spatial}$ and $L_{spectral}$ by Euclidean distance weights to make the central pixel contribute more to the reconstruction error of the patch. The two parts of the loss functions are described as

$$L_{spatial}={\displaystyle \sum }_{i=1}^p{\displaystyle \sum }_{j=1}^p{w}_{i,j}\ast \Vert x\left(i,j\right)-\widehat{x}{\left(i,j\right)\Vert }^2$$

(11)

$$L_{spectral}={\displaystyle \sum }_{i=1}^p{\displaystyle \sum }_{j=1}^p{w}_{i,j}\ast arccos\left(\frac{x\left(i,j\right)\widehat{x}\left(i,j\right)}{\Vert x{\left(i,j\right)\Vert }_2\Vert \widehat{x}{\left(i,j\right)\Vert }_2}\right)$$

(12)

where $x\left(i,j\right)$ denotes the spectral vector at the spatial location in row $i$ and column $j$ in the input cube of the original hyperspectral image, $\widehat{x}\left(i,j\right)$ denotes the spectral vector at the corresponding location in the reconstructed cube of that branch of the network, and $w_{i,j}$ denotes the weight of distance at that location, which is expressed as

$$w_{i,j}=\frac{1}{\sqrt{{\left(i-\lceil \frac{p}{2}\rceil \right)}^2+{\left(j-\lceil \frac{p}{2}\rceil \right)}^2}}$$

(13)

where $p$ is the spatial size of the 3D patch. Thus, the loss function of the first branch of the network is described as

$$L_1=L_{spatial}+\left(1-\alpha \right)L_{spectral}$$

(14)

where $\alpha$ denotes the weighting factor of the two parts; we set this as $\alpha =0.8$ in this study based on experiments with different values of $\alpha$. The loss function $L_2$ for the second branch of the network was calculated in the same way as $L_1$, except that $x\left(i,j\right)$ was replaced by $s\left(i,j\right)$. Thus, the loss function used for training the network can be described as

$$Loss=L_1+L_2$$

(15)

The training dataset is filtered using the background samples selection method described in Section 2.2.1 on the hyperspectral data introduced in Section 3.1. In the model training, the Adam optimizer [47] was applied. The training process continues until either 100 iterations have been completed or the loss function has converged. After training the network using the background samples, unfiltered hyperspectral data and multi-scale spectral difference data are used as input during the testing stage to obtain the reconstructed data. Since the network was only trained on background samples, the reconstruction error of anomaly patches should be very high. Therefore, model validation is performed by inputting unfiltered data into the network to generate a reconstruction error map. The results of spatial–spectral anomaly detection were obtained by calculating the spectral angular distance between the original hyperspectral image $X$ and its reconstructed data $\widehat{X}$, described as

$$A_{\mathbf{2}}=arccos\left(\frac{X\widehat{X}}{{\Vert X\Vert }_2{\Vert \widehat{X}\Vert }_2}\right)$$

(16)

We obtained the results of spatial detection $A_1$ and the results of spatial–spectral joint anomaly detection $A_2$, both of which measured the likelihood of an anomaly occurring in different dimensions. The final result of anomaly detection was obtained by fusing the two in the form of a product, which can be described as

$$A=A_{\mathbf{1}}\ast {\mathbf{A}}_{\mathbf{2}}$$

(17)

3. Experimental Setting and Results

In this section, the experiments used to evaluate the effectiveness of the proposed method are described in detail. First, we introduce the airborne and satellite-borne datasets used for the experiments. Then, the experimental setup is described, including the other methods in the comparison and the settings of the experimental parameters. Finally, we analyze and discuss the experimental results for the airborne and satellite-borne datasets.

3.1. Experimental Hyperspectral Datasets

Currently, most of the hyperspectral datasets used to validate anomaly detection methods contain airborne hyperspectral images, and fewer satellite-borne hyperspectral images are available. Therefore, in this study, we selected some images containing anomalous targets from a large number of hyperspectral images acquired by the GF-5 satellite AHSI to form a satellite-borne hyperspectral dataset named the G5 anomaly dataset.

In the experiments of this study, we selected four airborne hyperspectral images and four satellite-borne hyperspectral images from three datasets: the San Diego dataset, the Airport–Beach–Urban (ABU) dataset, and the G5 anomaly dataset. The details of the experimental images are as follows.

(1)

San Diego dataset: The San Diego dataset was acquired by the AVIRIS sensor. The first two airborne HSI used in the experiments are from this dataset. They have a spatial resolution of 3.5 m and a spectral resolution of 10 nm. The spatial size of the images is 100 × 100. After removing the water vapor bands and the low SNR bands, we selected 189 of the 224 bands with spectral coverage ranging from 370 to 2510 nm. The first image, denoted as SanDiego-I, has the anomalous target of three aircraft in the upper right corner, occupying a total of 58 pixels; the second image, denoted as SanDiego-II, has the anomalous target of three aircraft in the lower left and middle positions, occupying a total of 104 pixels. The pseudo-color images and ground truth maps of these two images are shown in Figure 5.

(2)

Airport–Beach–Urban (ABU) dataset: The third and fourth airborne HSI used in the experiment were from the ABU dataset. The third HSI was collected by the AVIRIS sensor. The spatial size of the image is 100 × 100, and the spatial resolution is 17.2 m. It has 198 bands selected from a total of 224 bands, with a range of 450–2500 nm and a spectral resolution of 10 nm. The anomalous target is the rock in the middle of the image in five columns, occupying a total of 155 pixels, denoted as Urban-I. The fourth HSI was collected by the ROSIS-03 sensor in Pavia, Italy. The spatial size of the image is 100 × 100, and the spatial resolution is 1.3 m. The number of bands is 102 ranging from 430–860 nm, with a spectral resolution of 3.3 nm. The anomalous targets are vehicles on the bridge, occupying a total of 68 pixels, denoted as Beach-I. The pseudo-color images and ground truth maps of these two images are shown in Figure 6.

(3)

G5 anomaly dataset: The AHSI on board the GF-5 satellite acquired a large number of valuable images [48], from which we selected images containing different anomalous targets in different scenes to establish a satellite-borne hyperspectral dataset for anomaly detection, named the G5 anomaly dataset. The images in this dataset are mainly from the visible near-infrared (VNIR) channel of the AHSI, with a spatial resolution of 30 m, a band number of 150, and a spectral resolution of 5 nm, with a spectral coverage ranging from 400 to 1000 nm. We selected four images of different typical anomalous targets in different typical scenes for the experiment. The images with a size of 100 × 100 pixels around the anomalous target were intercepted as the experimental images. G5-I is a pixel-level anomaly against a thin cloudy, and land background. G5-II is a building anomaly against a land background. G5-III is a ship anomaly against an ocean background. G5-IV is a building anomaly against a lake background. The pseudo-color images synthesized using the 74th, 38th, and 12th bands, and the ground truth maps are shown in Figure 7.

Figure 5. Airborne HSI from the San Diego dataset used in the experiment: (a) pseudo-color image of San Diego-I, (b) ground truth of San Diego-I, (c) false color image of San Diego-II, and (d) ground truth of San Diego-II.

Figure 5. Airborne HSI from the San Diego dataset used in the experiment: (a) pseudo-color image of San Diego-I, (b) ground truth of San Diego-I, (c) false color image of San Diego-II, and (d) ground truth of San Diego-II.

Figure 6. Airborne HSI from the ABU dataset used in the experiment: (a) pseudo-color image of Urban-I, (b) ground truth of Urban-I, (c) false color image of Beach-I, (d) and ground truth of Beach-I.

Figure 6. Airborne HSI from the ABU dataset used in the experiment: (a) pseudo-color image of Urban-I, (b) ground truth of Urban-I, (c) false color image of Beach-I, (d) and ground truth of Beach-I.

Figure 7. Satellite-borne HSI from the G5 anomaly dataset used in the experiment: (a) pseudo-color image of G5-I, (b) pseudo-color image of G5-II, (c) pseudo-color image of G5-III, (d) pseudo-color image of G5-IV, (e) ground truth of G5-I, (f) ground truth of G5-II, (g) ground truth of G5-III, and (h) ground truth of G5-IV.

Figure 7. Satellite-borne HSI from the G5 anomaly dataset used in the experiment: (a) pseudo-color image of G5-I, (b) pseudo-color image of G5-II, (c) pseudo-color image of G5-III, (d) pseudo-color image of G5-IV, (e) ground truth of G5-I, (f) ground truth of G5-II, (g) ground truth of G5-III, and (h) ground truth of G5-IV.

3.2. Comparison Algorithm and Parameter Settings

To evaluate the effectiveness of the method proposed in this study, six different state-of-the-art anomaly detection methods were selected for comparison. Among them, the statistics-based methods were RX [15] and LRX [49], the representation-based method was CRD [27], the decomposition-based methods were LRASR [31] and LSMAD [29], and the deep learning-based method was RGAE [35]. For each method, we set the parameters that optimized its performance. For LRX, the size of the sliding dual window ($W_{out}$, $W_{in}$) needed to be set. For CRD, the size of the sliding dual window ($W_{out}$, $W_{in}$) and the penalty factor $\lambda$ needed to be set. For LRASR, the number of clusters $\mathsf{{\rm K}}$, the number of pixels per class $P$, and the low-rank sparse trade-off factors $\beta$ and $\lambda$ needed to be set. For LSMAD, the maximum rank $r$ of the background matrix and the cardinality parameter $k$ of the sparse matrix needed to be set. For RGAE, the trade-off parameter $\lambda$ and the number of superpixels $S$ needed to be set. For our method, the dimensions of the hidden layer $d$ of the network, the morphological structure $se$, and the number of principal components $P$ needed to be set. For the four satellite-borne hyperspectral images from the G5 anomaly dataset, we set the same parameters. The parameters of the six methods and our method were set as shown in

Table 1. Main parameters for the different methods for different datasets.

Method	Parameter	San Diego-I	San Diego-II	Urban-I	Beach-I	G5
LRX	$(W_{out}$$, W_{in}$)	(25, 23)	(25, 23)	(21, 19)	(9, 7)	(9, 5)
CRD	$(W_{out}$$, W_{in}$)	(17, 15)	(17, 15)	(9, 5)	(9, 7)	(9, 5)
	$\lambda$	10−6	10−6	10−6	10−6	10−6
LRASR	𝐾	15	15	15	15	15
	𝑃	20	20	20	20	20
	𝛽	0.005	0.005	0.005	0.005	0.005
	𝜆	0.01	0.01	0.01	0.01	0.01
LSMAD	𝑟	2	2	2	1	2
	𝑘	0.005	0.005	0.01	0.01	0.002
RGAE	𝜆	0.01	0.01	0.01	0.01	0.01
	𝑆	150	150	100	150	150
Proposed	𝑃	2	2	2	2	2
	𝑠𝑒	(5, 5)	(5, 5)	(5, 5)	(5, 5)	(3, 3)
	d	200	200	200	200	200

.

The computer configuration used in the experiment had 256 GB of main memory, an Inter (IR) Xeon (R) Gold 5218 [email protected] GHz, and an NVIDIA GeForce RTX 3090 GPU, and the software resources used in the experiment include: Python 3.8.12 compiler, PyTorch 1.9.0 deep learning framework, and PyCharm IDE.

3.3. Experimental Results and Analysis

3.3.1. Experimental Results for the Airborne Hyperspectral Image Datasets

In this section, we describe and analyze the experimental results for the airborne hyperspectral images in detail. Figure 8 shows the results of detecting anomalies in the airborne datasets.

As can be seen in the results of visual detection (Figure 8b), our proposed method successfully highlighted the anomalous targets and suppressed the background well for all four datasets compared with the other six methods. The anomaly scores for the anomalous targets differed significantly from the background, and the anomalous targets could be clearly observed in the visual detection maps. This is because we made use of the spatial–spectral joint information along with the spatial information, which allowed for better suppression of the background and enhanced the response to the anomalous targets. For RX, the more complex the background of the HSI was, the greater the percentage of anomalous pixels and the worse the results of detection. As can be seen from Figure 8c, the RX algorithm failed to detect the anomalous targets and incorrectly detected many background pixels as anomalies in the San Diego-I dataset. For the other three datasets, the RX algorithm was able to detect some of the anomalies but still detected many background pixels as anomalies. For LRX, using pixels within a sliding dual window to estimate the background statistical features produced large errors and made it difficult to detect anomalous targets. Figure 8d shows that the LRX algorithm lost most of the anomalous targets in the San Diego-I, San Diego-II, and Urban-I datasets. For CRD, it can be seen from Figure 8e that the anomalous targets were detected in the San Diego-I, San Diego-II, and Beach-I datasets, but the difference between the anomalous targets and the background was small, with many background pixels responding more than the anomalous targets. The CRD algorithm lost some of the anomalous targets in the Beach-I dataset, probably because some of the lost anomalous targets were larger and denser.

For LRASR, Figure 8f shows that all anomalous targets could be clearly seen in all four datasets. However, its lack of suppression of the background caused the background to be too visible. LSMAD was able to separate the anomalies from the background very well. As can be seen in Figure 8g, the anomalous targets were clearly detected with less interference from the background and noise in the four datasets. For RGAE, as can be seen in Figure 8h, the results of detection were very close to those of our proposed method for the San Diego-II, Urban-I, and Beach-I datasets. However, for the San Diego-I dataset, the anomalous target was not clear enough, and the most obvious place was the background pixels.

Next, we analyzed the results of detection using three commonly used evaluation criteria: the receiver operating characteristic (ROC) curve, the area under the ROC curve (AUC), and the separability map. The ROC curve reflects the change in the FPR (false positive rate) and the TPR (true positive rate) for different thresholds, with the horizontal axis being the FPR and the vertical axis being the TPR. The closer the curve is to the top left corner, the better the result is. The AUC value is the area under the ROC curve, and the closer it is to one, the better the result is. The separability map shows the dispersion of background and abnormal pixels in the results and reflects the ability of the results to separate the background and the anomalies.

Figure 9 shows the ROC curves of the different methods for the four airborne HSI. We can see that in Figure 9a, the curve of our proposed method was always at the top for the San Diego-I dataset. This means that our method had the maximum detection rate for different false positive rates. For the San Diego-II dataset in Figure 9b, the curve of our proposed method had a high detection rate at lower false positive rates, while the RGAE method was able to achieve the maximum detection rate first. As can be seen in Figure 9c,d, for the latter two datasets, our proposed method was the first to achieve a high detection rate at a lower false positive rate.

Table 2. AUC values of the results of detection of different methods for each airborne HSI.

AUC	San Diego-I	San Diego-II	Beach-I	Urban-I
RX	0.9053	0.9403	0.9885	0.9951
LRX	0.8725	0.9675	0.9284	0.9188
CRD	0.9788	0.9293	0.9570	0.9283
LRASR	0.9836	0.8803	0.9778	0.9456
LSMAD	0.9864	0.9813	0.9903	0.9927
RGAE	0.9791	0.9919	0.9914	0.9887
Proposed	0.9974	0.9927	0.9940	0.9980

shows the AUC values of the results of different detection methods for each airborne HSI. It can be seen that for airborne datasets, our proposed method had the largest AUC value compared with the other six methods.

Figure 10 shows the separability map of the different methods for the four airborne HSI. It can be observed that our proposed method was able to achieve an effective separation of the anomalies and the background in the four airborne HSI. The proposed method was able to achieve high anomaly scores while also effectively suppressing the background. The LRASR and RGAE methods were equally capable of separating anomalies from the background in the four datasets, but neither suppressed the background as well as our proposed method.

3.3.2. Experimental Results for the Satellite-Borne Hyperspectral Image Dataset

We carried out the same experiment with the four satellite-borne HSI. Figure 11 shows the visual detection results. It can be seen that the RX method was able to highlight most of the anomalous targets, while in the case of complex background components, it highlighted a large amount of background information at the same time, such as in G5-II, G5-III, and G5-IV. Compared with the RX method, the LRX method suppressed the background better, while the reinforcement of the anomalous target became weaker. For anomalous targets with a large spatial area, such as G4, the detection was poor. The RX and LRX method achieved better results with G5-I because the distribution of the image background was more uniform and was close to a Gaussian distribution. The CRD method was weak at suppressing the background, and the background pixels generally had a high anomaly score. The LRASR method highlighted the anomalous targets very well in all four HSI; however, the background was not well suppressed in all four HSI. For the LSMAD method, when the background had a small number of components and a homogeneous distribution, such as in G5-I, G5-II, and G5-III, this method was able to detect the anomalous targets and suppress the background. For G-IV, which had more complex backgrounds, it lost some of the anomalous targets and had high scores for more background pixels. The RGAE method was able to highlight most of the anomalous targets in the first three HSI, but the background was poorly suppressed; in the last image, it lost some of the anomalous targets. Our proposed method clearly had the best ability to highlight the anomalous targets in all four HSI. Due to the combination of spatial and spectral information, our proposed method was able to suppress the background better; only a few high-contrast edge background pixels had high anomaly scores.

The ROC curves and AUC values for the different methods for the four satellite-borne HSI are presented in Figure 12 and

Table 3. AUC values of different methods for four satellite-borne HSI.

AUC	G5-I	G5-II	G5-III	G5-IV
RX	0.9999	0.9945	0.9135	0.9663
LRX	0.9999	0.9903	0.8865	0.8691
CRD	0.9999	0.9938	0.8758	0.6394
LRASR	0.9999	0.9992	0.9993	0.9720
LSMAD	0.9997	0.9995	0.9998	0.5482
RGAE	0.9997	0.9619	0.9929	0.6136
Proposed	0.9999	0.9995	0.9998	0.9949

. It is important to specify that for G5-I, where the anomalous target occupied a few pixels, RX, LRX, CRD, LRASR, and our proposed method all obtained the highest scores for the anomalous target’s location. Thus, as shown in Figure 12a, the ROC curves were all red horizontal straight lines with a value of one. These methods had an AUC value of 0.999. For the other three HSI, the ROC curves of our proposed method were higher than those of the other methods and achieved a detection rate close to one with the lowest false positive rate. The proposed method achieved the highest AUC values for these four HSI.

For the satellite-borne HSI with anomalous pixels occupying only a few pixels, the ROC curve was close to a horizontal straight line with a value of one, and the AUC value was close to one when the anomalous pixels had the highest score, such as G5-I. In this case, the ROC curve and AUC values only indicated that anomalies can be detected but did not reflect the effect of suppressing the background. Therefore, as shown in Figure 13 and Figure 14, we used a 3D mesh map and a separability map to visualize the suppression of the background and the degree of separation of the anomalous targets and the background. It can be seen that for G5-I, although the other methods showed good performance in terms of the ROC curves and AUC values, their ability to suppress the background was weaker than that of our proposed method. For all four satellite-borne HSI, our proposed method adequately suppressed the background while separating the background from the anomalous targets.

Based on the detection results of airborne and satellite-borne hyperspectral images shown in Figure 8 and Figure 11, we discuss the differences in anomaly detection results between airborne and satellite-borne hyperspectral images. Satellite-borne hyperspectral images have a larger spatial range compared to airborne hyperspectral images, but the spatial resolution is somewhat lower. Therefore, for some small anomalous targets, satellite-borne hyperspectral images may exhibit pixel-level anomalies, which require anomaly detection methods to make full use of spectral information, such as Figure 11a. In addition, although satellite-borne hyperspectral images have lower spatial resolution due to their large spatial range, the background distribution in a local area may be more uniform compared to airborne hyperspectral images, resulting in better anomaly detection results. Our proposed method can effectively utilize both spatial and spectral information of hyperspectral images and achieve good detection results on both airborne and satellite-borne hyperspectral images.

4. Conclusions

In this study, we proposed a novel method of detecting hyperspectral anomalies based on the two-branch 3D convolutional autoencoder and spatial filtering with the aim of making full use of the spatial–spectral joint information. The proposed method consists of two parts: spatial detection and spatial–spectral joint detection. During spatial detection, we extracted the spatial features using a morphological filter and a total variational curvature filter. In spatial–spectral joint detection, a two-branch 3D convolutional autoencoder network was proposed to extract the spectral interband features and spatial–spectral joint features from hyperspectral images. Due to the lack of satellite-borne hyperspectral images for evaluating anomaly detection, we constructed a satellite-borne hyperspectral dataset containing different backgrounds and different anomalous targets based on the HSI acquired by the GF-5 satellite. To verify the effectiveness of our proposed method, experiments were conducted on four airborne hyperspectral images and four satellite-borne hyperspectral images with different backgrounds. The experimental results showed that compared with other state-of-the-art methods, our proposed method performed well on both airborne and satellite-borne datasets, highlighting anomalies while suppressing the background. However, we have not yet conducted experiments on whether our two-branch 3D-CAE model can be transferred to other fields, such as low-rank recovery, classification, and denoising. We will conduct further research in these areas in the future.