Hyperspectral Anomaly Detection With Morphological Random Walker

This work proposes a new morphological random walker (MRW) method for hyperspectral anomaly detection. The proposed method introduces a morphology-based objective function into a random walker (RW) algorithm, sufficiently exploiting spatial morphological property and spatial similarity of HSIs for detection. The MRW method comprises two major stages. Firstly, we employ the extended morphological profiles (EMPs) and different operations to extract the spatial morphological property of HSIs. Second, according to the morphological property, we construct a morphology-based objective function. This function is incorporated into the RW-based optimization model, encoding the spatial similarity of HSIs in a weighted graph. Two factors determine the class of test pixels, including the spatial morphological information learned by EMPs, and the spatial correlation among adjoining pixels modeled by the weighted graph. Since the two factors are well considered in the MRW method, the proposed method illustrates outstanding detection performances for several widely used real HSIs.

Over the past decades, numerous methods have been developed for anomaly detection in HSIs [15]- [18]. In general, we can divide these detection approaches into three major categories. One category of approaches is based on statistical modeling [19]. Statistical modeling techniques have a close relationship with the Gaussian distribution. The well-known Reed-Xiaoli (RX) method [19] employs a The associate editor coordinating the review of this manuscript and approving it for publication was Shuihua Wang . multivariate Gaussian model to characterize the background information. Anomaly targets are detected by estimating the Mahalanobis distance of a pixel vector to the background. Furthermore, several improved RX-based methods have been proposed [20]- [25], such as weighted-RX [21], kernel-RX [22], and cluster kernel-RX (CKRX) [24]. These methods usually can generate outstanding detection results. Moreover, several robust statistical modeling techniques have been proposed, including robust no-linear learning-based detection method [26], and random-selection-based detection method [27].
Instead of estimating statistical models, many researchers pay attention to geometrical modeling methods [30] for hyperspectral anomaly detection. These methods judge that background pixels can be constructed by some main pixels, whereas anomalies cannot. The sparse representation-based [31]- [33], collaborative representation-based [34], [35], and the low-rank representation-based detectors [36]- [38] have been developed in recent years. The collaborative representation detector (CRD) [34] captures the difference between a test pixel and its surrounding pixels. The background joint sparse representation detection (BJSRD) method [31] is proposed to utilize the redundant background information and deal with the complicated multiple background classes. Moreover, a spectral unmixing and dictionary-based low-rank decomposition method [37] explores the low-rank prior information of background, and distinguishes the sparse anomalies from the low-rank background. Moreover, some deep learning-based detectors [39]- [41] are also employed to detect anomalies. These methods can separate targets from the background, by extracting high-level features. The detection results are generated by constructing a trained convolution neural network (CNN) model [39].
These methods have achieved an effective detection performance by exploring spectral-spatial information in HSIs. Therefore, in this paper, the motivation is that how to make more efficient use of spectral-spatial information in HSIs for hyperspectral anomaly detection. In the case of scene classification, numerous spectral-spatial feature extraction methods [1], [44] have been proposed to enhance the classification performance via capturing morphological property in images, such as attribute profiles (AP) [45], [46], extinction profiles (EP) [47], [48], morphological profiles (MP)based methods [49], [50], and its extended versions (EMP) [51], [53]. Among these methods, the EMP has drawn lots of attention. Because it could fully exploit the spatial morphology property. In addition, the random walker (RW) algorithm [54], [55] has been widely used in hyperspectral image classification, since it can effectively explore spatial similarity information of the image. For example, Kang et al. [56], [57], and Cui et al. [58] develop RW-based methods to explore the spectral-spatial information of HSIs for classification. The RX-based methods usually require the user to label some target samples. It means they may not be employed directly in the unsupervised anomaly detection field.
To exploit spectral-spatial information in anomalies effectively, this paper introduces a morphological random walker (MRW) technique, jointly capturing the spatial morphological property and spatial similarity information in HSIs for hyperspectral anomaly detection. The MRW technique comprises two major stages. First, EMPs and different operations extract the spatial morphological property in HSIs. Then, according to the morphological property, we construct a morphology-based objective function, and incorporate this function into the RW algorithm for detection. The main contribution of this paper is to propose an effective anomaly detection method, jointly capturing both the morphology property of anomaly objects and the spatial similarity among adjoining pixels simultaneously.
The rest of this work is presented as follows. We review the EMP feature extraction method and the RW algorithm in section II. We introduce the MRW method in Section III. Experimental results are presented in Section IV. Last, conclusions are given in Section V.

II. RELATED WORKS A. EMP
A morphological profile (MP) [49] is constructed by employing a series of morphological opening and closing with structuring elements (SEs) on the principal components (PCs). When MPs on several PCs are integrated, the EMP can be generated. Specifically, there are two fundamental morphology operators, i.e., erosion and dilation. Opening and closing are combinations of erosion and dilation. However, two operations may change some intrinsical structures in the image, which may produce some fake objects. In [51], [52], geodesic morphology and reconstruction are introduced, which make opening and closing satisfy the following assumption. If the structure of the image fails to contain the SE, it would be removed. Conversely, it would be preserved. Therefore, it is necessary to employ SEs in different sizes. A MP consists of the opening profile (OP) and the closing profile (CP). We define the OP at the pixel x of the image I as a n-dimensional vector, where γ (i) R means the opening by the reconstruction with the SE of the size i, and n is the total number of openings. Similarly, the CP at the pixel x of the image I is also a n-dimensional vector, where φ Therefore, an image can generate a multiband image. We note that the MP is just constructed by employing a PC in HSIs, which may lose some spectral information. In [51], it is suggested to use several PCs of the HSI, which contains the most important spectral information in the HSI. When MPs on several PCs are integrated, the EMP feature will be generated. The EMP is a [m(2n + 1)]-dimensional vector, where m is the number of retaining PCs.

B. RW
The RW algorithm is originally designed for image segmentation [54], [55]. This algorithm can group the pixels of the image X into various classes [56]. First, some pixels Q are labeled from X . These pixels serve as the seeds, with at least one pixel of each class. The seeds are defined as , where x T M and x T U are the marked seed nodes and the unmarked nodes, respectively).
Then, the RW model defines an image as a graph G = (V , E) with vertices v ∈ V and edges e ∈ E. Vertices are pixels in the image, and edges contact two adjoining pixels. Each edge e ij links the ith and jth pixels based on a weight ω(e ij ) = e −β(v i −v j ) 2 , and β is a parameter. We define a n × n Laplacian matrix L as, , if i and j are adjacent pixels. 0, otherwise Then, let p l = [p l 1 , . . . , p l n ] represent the probability vector of X for class l. It is decomposed as represents the seeds having a fixed value 0 or 1. p l can be obtained by minimizing the Dirichlet integral [55]: We differentiate Dir[p l ] with p l U , and the critical point is found as: With the labeled pixels p l M , we can calculate the energy function shown in (6), by solving a series of linear equations [54]. Fig.1 displays the flowchart of the MRW method. It consists of two major stages. First, EMP and differential operations capture the spatial morphological property of HSIs. Then, according to the spatial morphological property, we construct a morphology-based objective function and incorporated it into the RW algorithm, encoding the spatial similarity of HSIs in a weighted graph. The detection result is obtained by jointly exploiting the spatial morphological information learned by EMPs and the spatial correlation among adjoining pixels modeled by the weighted graph.

A. EXTRACTION OF THE MORPHOLOGICAL PROPERTY OF ANOMALIES
First, the spectral dimension of an input HSI I is reduced by performing the principal component analysis (PCA) [61]. The EMP of the I is obtained as follows, where m is the number of retained PCs. An example of MP is illustrated in Fig. 1 (the number of opening/closing is set as 1 and the corresponding structure size setting to 2). Given an HSI [see Fig. 2(a)], we estimate the MP of its first PC. The first PC of HSI and the resulting morphological operations are shown in Fig. 2(b)-(d). We use the opening and closing transforms to separate bright (see Fig. 2(c)) and dark (see Fig. 2(d)) structures from the image. The bright/dark is brighter/darker than the adjoining pixels features [52]. Then, to capture the spatial morphology property of anomalies, the differential operations compute the differences between each MP and the corresponding PC [60]. In Fig. 2(e) and (f), due to the difference between anomalies and background pixels on bright and dark structures, the morphology property of anomalies (such as shape, size, and location) could be explored. The differential map s emp is a 2nm dimensional vector, where s mp PC s represents the sth differential map in the sth PC. |OP n(PC s ) − PC s | is the morphological feature of the bright objects (see Fig. 2(e)), and |CP n(PC s ) − PC s | is the morphological feature of the dark objects (see Fig. 2(e)). To obtain the fused morphology features P = (P 1 , P 2 ), the differential maps are fused, in which P 1 i ∈ [0, 1] represents the probability that the ith pixel pertains to the anomalies.
where P 1 represents the initial detection map (e.g., in Fig. 2(g)), effectively reflecting the spatial morphology property in HSIs.

B. MRW-BASED OPTIMIZATION
To further capture spectral-spatial information in HSIs for hyperspectral anomaly detection, we construct a morphology-based objective function, and incorporate this function into the RW model. Thus, in this letter, the proposed algorithm could fully exploit the spatial morphological property in HSIs and spatial similarity information of adjoining pixels. Specifically, we construct a new morphology fitting , restricting the Dirichlet integral [54] to be as close to the prior anomaly distribution as possible [60]. Then, we incorporate the morphology fitting constraint into (4), and the final detection result is estimated as follows, where µ denotes a controlling parameter. The second term is a morphology fitting constraint. Y means a pixel-wise indication vector, which inherits the values of P. The proposed method is estimated in a pixel-wise manner. p l and Y are v×1 vectors. L means a v × v matrix. v is the number of pixels in the image.
To generate the first term (spatial similarity term), we construct a weighted graph G = (V , E) based on the first PC of the HSI [60]. In this graph, V represents pixels in the first PC. E means edges that link the pixels. ω(e ij ) = e −β(v i −v j ) 2 is defined for each edge e ij to characterize the intensity difference between the adjoining pixels in G. To generate the second term (spatial morphological term), the initial detection map P is employed to estimate the pixel-wise indication vector in (12). The morphological fitting constraint in (12) can provide a prior morphology property of HSIs. Therefore, it could better generate detection performance. To choose seeds from the P 1 , a threshold t high is defined, which are utilized to select pixels with Y U > t high as anomalies seeds. These seeds are defined as p l M . Thus, the matrix decomposition of (12) is as follows, Since the seeds are generated from the initial detection map automatically, MRW needn't user interaction, which has more practicability in hyperspectral anomaly detection.
Then, p l M and p l U are combined as p l . p 1 represents anomaly possibility, and we reshape it to a matrix S final . The same size of the input image as the final detection result.  In Fig. 2(h), errors in the anomalies and background areas can be removed effectively. Moreover, the boundaries between anomalies and background areas become more accurate. The reason is that, aside from the morphological property in HSIs, the spatial correlation among adjacent pixels is also considered in the MRW-based method.

IV. EXPERIMENTS AND DISCUSSION
A. DATA SETS Some real hyperspectral data sets, including the San Diego, ABU-Beach, and ABU-Urban data sets [18], are employed to assess the performance of the MRW detector. These images can be download on this website. 1

1) SAN DIEGO DATA SET
This data set is captured by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) sensor. It is over the San Diego airport area, CA, USA. This data set comprises 100 × 100 pixels with a spatial resolution of 3.5 m per pixel. It has 189 spectral channels in the wavelength range of 0.37-2.51 um. In this data set, we consider three airplanes as the anomaly targets. The test image and the corresponding reference detection map are displayed in Fig. 3.

2) ABU-BEACH AND ABU-URBAN DATA SETS
The second and third data sets include three test images respectively. These images are manually extracted from the   AVIRIS web site. 2 In Table 1  To estimate the effectiveness of the proposed MRW detector, five widely used detection methods, RX method [19], local RX (LRX) method [19], local kernel-RX (LKRX) method [22], LRaSMD-based Mahalanobis-distance anomaly detection (LSMAD) method [36], and attribute and edge-preserving filtering-based anomaly detection (AED) method [18] are used to compare. The detection results are accessed by the receiver operating characteristic (ROC) area under the curve (AUC) metric [33]. ROC curves are generated based on true positive rate (TPR) and false positive rate (FPR) with a threshold range of [0, 255]. At each possible threshold, the detection result is binarized into 1 and 0 to mean the anomaly target and background. A better detector usually lies nearer the upper leftmost corner [63]. Moreover, we also computed the AUC score to assess the performance of these detectors. Generally speaking, an outstanding detection model can obtain a high AUC value [64]. 2 http://ariris.jpl.nasa.gov

B. PARAMETERS DISCUSSION
For the MRW detector, the main parameters include the number of PCs, i.e., m, and the parameters of the optimization, i.e., β and µ. The AUC score is used to evaluate the detection performances of the proposed method under different parameter settings. To generate the EMP, four openings and closings are estimated for each PC, and the corresponding structure sizes are set to 2, 4, 6, and 8. Fig. 6 shows the AUC value of the MRW detector with different values of β and µ. It can be observed that µ should be lower than 10 −3 and higher than 10 −5 , which can generate the outstanding performance. Moreover, β should be higher than 30. It is obvious that ω(e ij ) = e −β(v i −v j ) 2 for adjacent pixels tend to be 1 if β is small. Thus, the edge weights fail to model the spatial similarity among adjoining pixels. Conversely, ω(e ij ) will be close to 0 when β is large. Thus, µ should be small to ensure the balance of two terms in (12). Otherwise, the morphology term will play a significant role in the detection, if µ is very large. In this situation, the final detection result would be similar to the initial detection result. In contrast, the spatial effect may cause the over-smoothed detection map if µ is too small. Therefore, β = 50 and µ = 10 −3 are set as the default parameters for the experiments.   Moreover, the number of PCs also is a significate parameter for the MRW detector. Fig. 7 illustrates the average detection scores of the MRW method on three data sets, i.e., the San Diego, ABU-Beach, and ABU-Urban data sets. As can be observed, when PCs = 1, the accuracies of the proposed method are relatively low. The useful spectral information will be lost if the number of features PCs is small. In contrast, the detection performance will decrease if PCs is larger than 5. This is because redundancy information will be produced if the number of features PCs is large. We can find that 2 ≤ PCs ≤ 4 can obtain satisfactory AUC value for all images. Therefore, we set PCs = 3 as the default parameter.

C. DETECTION PERFORMANCE
Experiments are performed on the San Diego, ABU-Beach, and AUB-Urban data sets. Figs. 8-10 show the detection maps obtained by different detection methods. For the LRX and LKRX algorithm, the window size w in is set to [3,19], and w out is [5,13]. For the LSMAD algorithm, the cardinality k is fixed to 0.005, and the rank r is changing from 1 to 20. Furthermore, for the AED algorithm, two filtering parameters σ s and σ r are set to 5 and 0.5 respectively.
For the San Diego data set, the obtained detection results are displayed in Fig. 8. The proposed MRW method gives a map where the anomalies are apparent. The RX method fails to detect anomalies effectively. The LRX detector has a better detection result than the RX detector. The LKRX method can well detect the positions of anomalies, but the shape of the objects is missing. For the LSMAD method, a part of background pixels are considered as anomaly targets. The MRW and AED methods both can achieve good detection results. The MRW method obtains more accurate boundaries between anomalies and background regions, since the AED method is sensitive to parameters for image smoothness. Therefore, some boundary regions may be over-smoothed when the parameters are inaccurate. Besides, the AUC scores are provided in Table 3. The proposed MRW detector achieves a score that is 5.2% higher than the RX detector, which is a distinct improvement. It also confirms that the proposed method can outperform the traditional detectors.
A main advantage of the MRW method is that its outstanding detection performance in different scenes. For the ABU-Beach image, the obtained detection results are shown in Fig. 9(a)-(c). The MRM model can provide a distinguishable detection map, highlighting the anomalies in different scenes. The AUC scores are illustrated in Table 4. It can be concluded that MRW is a promising method for anomaly detection in complex scenes.
The MRW method also can achieve outstanding detection performance, when HSIs contain serious noise. For the ABU-Urban data sets, which are influenced by the serious strip noise (see Fig. 5), the obtained detection maps are illustrated in Fig. 10(a)-(c). In Fig. 10(a), it is obvious that the RX and LSMAD methods obtain a better detection performance than LRX and LKRX methods. Specifically, the LRX method losses some anomaly pixels. Meanwhile, the LKRX method fails to detect anomaly targets effectively. The corresponding AUC scores are provided in Table 5. Although the AED method can achieve higher detection scores than the RX and LSMAD methods, it fails to reduce noise and suppress its disturbing effectively. Compared with the AED method, it is evident that the MRW method can perform better at noise pixel suppression. The reason is that the proposed MRW method can effectively exploit spatial correlations between adjoining pixels to reduce the interference of noise.
Then. we further discuss the influence of the MRW-optimization step in the proposed method. Specifically, the LRX detector is adopted to processing the initial  [19], (b) LRX method [19], (c) LKRX method [22], (d) LSMAD method [36], (e) AED method [18], (f) MRW method.   detection result in Section III-A. Fig. 12 displays the detection maps generated by the DEMPs, LRX, MRW, and DEMPs + LRX methods, respectively. DEMPs represents the initial detection result in the proposed method. It can be observed that the MRW and DEMPs + LRX methods produce better detection performance than others without   relying on the spatial morphology information in HSIs. Fig. 13 illustrates that the MRW method has AUC = 0.9908, which is better than others.
The corresponding ROC curves of different methods are shown in Fig. 11. We can observe that the MRW method usually has higher true positive rates. The computational complexity of the proposed method is as followed. The main running time is consumed by constructing the EMP feature, performing differential operations, and estimating the MRW optimization. Specifically, the EMP feature is constructed by integrating MPs on several PCs. Each MP feature has (2n + 1) band images. Thus, the computational complexity of constructing the EMP feature is o ((2n + 1)m). m is the number of retaining PCs. Then, the computational complexity of differential operations is o((2n + 1)σ ). σ means the differential operation for each band in the EMPs features. Last, we denote the number of pixels in HSIs is v. For the MRW optimization step, the computational complexity is based on the total input pixels v, and the time-complexity of evaluating anomaly scores is o(vτ ). τ represents the Dirichlet integral operation in Eq. (14). Therefore, the computational complexity of the proposed method is  o((2n + 1)(m + σ ) + vτ ). Furthermore, the running time of different methods is reported in Table 6 (For AUB-Beach and AUB-Urban data sets, the computing time is based on calculating the average time of the corresponding three subgraphs). The programs of various methods are executed on a computer with a 2.6-GHz CPU and 32-G memory, and the software platform is MATLAB R2017b. From Table 6, we can observe that the LRX and LKRX models take more running time. Although the MRW method is not the most efficient among the compared ones, it can obtain much better detection results.

V. CONCLUSION
In this work, a novel MRW detector is presented for hyperspectral anomaly detection. The MRW detector consists of two steps, i.e., extraction of the morphology property of anomalies and MRW-based optimization. First, the extended morphological profiles and differential operations are applied to extract the spatial morphological property of HSIs. Then, according to the morphological property, we construct a morphology-based objective function and incorporate it into the RW model for detection, jointly exploiting the morphological property of anomalies and the spatial correlations among adjacent pixels. Experimental results demonstrated that the MRW algorithm could achieve outstanding detection performances on three real HSI data sets.
The parameters empirically setting way is a limitation of the proposed detector when it is employed in reality detection applications. Therefore, automatically choose optimal parameters is the focus of our future works. Besides, applying MRW in other HSI applications will also be further investigated (i.e., segmentation, visualization, and change detection).