Hyperspectral anomaly detection based on the bilateral filter

Accepted Manuscript Hyperspectral Anomaly Detection Based on The Bilateral Filter Xifeng Yao, Chunhui Zhao PII: DOI: Reference:

S1350-4495(18)30186-5 https://doi.org/10.1016/j.infrared.2018.05.028 INFPHY 2579

To appear in:

Infrared Physics & Technology

Received Date: Revised Date: Accepted Date:

21 March 2018 26 May 2018 30 May 2018

Please cite this article as: X. Yao, C. Zhao, Hyperspectral Anomaly Detection Based on The Bilateral Filter, Infrared Physics & Technology (2018), doi: https://doi.org/10.1016/j.infrared.2018.05.028

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Hyperspectral Anomaly Detection Based on The Bilateral Filter Xifeng Yao and Chunhui Zhao

 Abstract—Bilateral filter outputs a weighted average of the neighboring information at a pixel. It can smooth the image and remove anomalous points while preserving edges. Inspired by this idea, this paper proposes a hyperspectral anomaly detection algorithm based on the bilateral filter. A pixel that spectrally differs from the background is regarded as anomalous. If every pixel in a hyperspectral image is represented by its surrounding background pixels, a relatively clean hyperspectral background image without anomalies can be obtained. We improve the weights that are imposed on the neighboring pixels in representation. First, spatial distance weight is selected to make full use of the space-varying information. Second, spectral Euclidean distance weight is used to express the spectral similarity. To reduce the negative impact that is caused by anomalies existing in neighboring background pixels, dual window is utilized to act on above two weights. Finally, anomalous information of hyperspectral image can be derived after making the subtraction between the original hyperspectral image and the derived hyperspectral background image. Extensive experiments on three hyperspectral images demonstrate that the proposed method outperforms other benchmark competitors. Index Terms—Hyperspectral remote sensing, anomaly detection, bilateral filter, weighted average, clean background representation.

I. INTRODUCTION imagery (HSI) processing has been received a wide and increasing attention for many years. The captured H yperspectral hyperspectral data are a cube, and it provides almost continuous spectral curves of ground materials on its third dimension [1], [2]. With these hundreds and even thousands of spectral bands, hyperspectral remote sensing can possess an intrinsic advantage to find and identify the subtle differences of different objects on the ground surface [3], [4]. Anomaly detection (AD), as a very important and fundamental branch of HSI, has been popular for many years [5]. In essence, AD is a binary classification problem, which is aimed to distinguish an anomalous point from the background at a pixel [6]. Different from the supervised target detection, hyperspectral AD is an unsupervised technique and does not require any prior information about target or background. Thus, this technique is more suited to the practical situation, in which the spectral information of many ground materials cannot be known as priors in advance [5], [7-10]. In this case, AD can be applied to detecting many different kinds of anomalous targets, such as special species in agriculture and ecology, rare minerals in geology, toxic wastes in environmental monitoring, oil spills in water pollution, drug trafficking and smuggling in law enforcement, man-made objects in battlefields, unusual terrorism activities in intelligence gatherings and tumors in medical imaging [11], [12]. The current hyperspectral anomaly detection algorithms are mainly to find discriminative differences between the pixel under test and the background, and then they estimate anomaly degree of the testing pixel by virtue of these differences [13-15]. So how to accurately treat and take advantage of the background has a very significant and vital impact on achieving a desirable detection performance. According to the processing manner of background, the existing anomaly detection algorithms can be divided into two categories. One category assumes that the hyperspectral background conforms to a Gaussian normal distribution, while the other one does not require such an assumption [16], [17]. The well-known RX detector (RXD) [18] proposed by Reed and Yu is a typical type of the first kind, and it is usually deemed to be the benchmark anomaly detection method among all [19]. RXD builds the two-classification hypothesis problem and designs the decision function from a generalized likelihood ratio test for an unknown additive contrast signal in a multivariate Manuscript received —; revised —; accepted. Date of publication — ; date of current version — . This work was supported by the National Natural Science Foundation of China (No. 616011351 and No. 61571145), and the Fundamental Research Funds for the Central Universities (No. HEUCF1608). Xifeng Yao was with the College of Information and Telecommunications, Harbin Engineering University, Harbin, 150001 China. He is now with the Space Science and Engineering Center, University of Wisconsin-Madison, Madison, USA (e-mail: [email protected]). Chunhui Zhao is with the College of Information and Telecommunications, Harbin Engineering University, Harbin, 150001 China (e-mail: [email protected]).

1

Gaussian background. RXD usually employs a sliding window to construct the local RX detector (LRXD), which has a capability to detect the local anomalies that are overwhelmed in the global RX detector (GRXD) [19], [20]. However, both LRXD and GRXD compute the background statistics by making use of all pixels in the local region or the whole image. Thus, the anomalies involved in the background can damage the Gaussian normal distribution, which reduces the separability between anomalous targets and background [21], [22]. To overcome this issue, many researchers have made efforts and contributions. Billor et al. [23] proposed a blocked adaptive computationally efficient outlier nominator (BACON) anomaly detector. This method tries to select the growing subset of outlier-free points and obtain robust background statistics to weaken the contamination by anomalies. The random-selection-based anomaly detector (RSAD) [21] proposed by Du and Zhang utilizes a random selection process to lower the sensitivity of contamination. The probabilistic anomaly detector (PAD) [24] obtains potential anomalous targets and background clutter based on the results of RXD. The median-mean line RX detector (MML-RXD) [25] is employed to provide more reliable background samples. However, it is hard for the hyperspectral background to satisfy the Gaussian normal distribution in reality, even though a background purification method is used [13], [26]. In addition, most of the time the algorithms designed by the Gaussian normal assumption have to perform the inverse operation of background covariance matrices to “whiten” the background statistics [27]. And these covariance matrices are generally singular, so some performance can be diminished [28], [29]. For the second type of AD methods without the Gaussian normal assumption, they usually do not have the restriction that the hyperspectral background follows a Gaussian normal distribution. Low rank and sparse matrix decomposition (LRSMD) [30] presented by Sun et al. is a good example. It is considered that hyperspectral background is low-rank and anomalies are sparse [31], [32]. LRSMD decomposes the whole hyperspectral image into a low-rank background matrix and a sparse anomaly matrix via the singular value decomposition method. And then it captures the anomalous information from the sparse matrix. The collaborative-representation-based detector (CRD) [33] makes an endeavor towards the conception that a background pixel can be represented by its surrounding neighbors while an anomalous pixel cannot. This method uses local background pixels to linearly fit the testing pixel by l2 -norm minimization of the representative weight vector and accomplish the anomaly detection through a distance function. Unfortunately, the inverse operation of matrices is also existed in CRD. In this paper, we propose a bilateral-filter-based anomaly detector (BFAD). Interestingly, BFAD belongs to the second anomaly detection algorithms without a Gaussian normal assumption, and it does not need a requirement that the matrices have to be inverted. This method is inspired by the bilateral filter (BF) in digital image processing [34], [35]. Bilateral filter substitutes the weighted surrounding neighbors for centric pixels to remove abnormal points. It should be noted that the processing objects of bilateral filter in digital image processing are gray images or color images. For the three-dimensional hyperspectral cubes, we should make full use of the information that they provide to highlight the difference between anomalies and background. So two types of weights including spatial distance weight (SDW) and spectral Euclidean distance weight (SEDW), are selected to construct the hyperspectral background image without anomalies. The dual window structure is a necessity in construction, which can increase the probability that the neighboring information used for constructing background pixels is pure. Finally, anomalous information can be found in the difference between the original hyperspectral image and the derived hyperspectral background image. The structure of the paper is as follows. Section II detailedly describes the proposed algorithm of BFAD. Next, in section III, three hyperspectral datasets are conducted to prove the efficiency of the proposed algorithm. Finally, the conclusions are summarized in section IV. II. ANOMALY DETECTOR BASED ON THE BILATERAL FILTER Bilateral filtering is very typical and popular in digital image processing, which is owing to the characteristic that it can smooth images while preserving edges [34]. Generally, image filtering is the process of a weighted average, and each pixel in the filtered image is derived from an average of its weighted neighboring pixels in the original image. Bilateral filter considers both geometric closeness and photometric similarity between pixels when choosing the reasonable weights. Thus, it can smooth colors and preserve edges in a way that is tuned to human perception [36]. And a weighted average has another advantage that the noise corrupting the neighboring pixels can be mutually less correlated and thus averaged away. In this work, we make a try to apply the idea of bilateral filter to hyperspectral anomaly detection. What it differs from gray images and color images lies in that every pixel in a hyperspectral image denotes a spectrum of the ground material(s). Thus, it is rational that hyperspectral background image can be obtained via the thought of a weighted average of the surrounding pixels. Finally, an anomaly decision function is operated on both the original hyperspectral image and the derived hyperspectral background image to accomplish the anomaly detection. A. Spatial Distance Weight In digital image processing, Gaussian low-pass filtering computes a weighted average of neighboring pixel values, in which the weights are only dependent on the distance from the neighborhood center and take on a downward trend [34]. It is generally 2

considered that images are typically varying slowly over space, so a pixel can be represented by its neighboring pixels. And this slowly varying forms an idea that the pixels which are closer to the center pixel will have more similar related features, while there will be more differences occurred in those pixels which are far from the center pixel. The above conclusions are also applicable to hyperspectral remote sensing images. In [37], Zhao et al. believe that the background spectral signal is stable and varying slowly. And according to this slowly varying feature, four kinds of anomaly estimators are presented to distinguish the anomalies. Therefore, we make full use of the space-varying information and give the first kind of weight, spatial distance weight (SDW). Let x and  denote the intensity of center pixel (i, j ) and any pixel

(k , l ) in the neighborhood, then the weight imposed on  can be expressed by:  1  d ( , x )  2  SDW (, x)  exp    d    2  d    

(1)

d d (, x)  (i  k ) 2 +(j  l ) 2 where dd (, x) measures the geometric closeness between the neighborhood center x and a nearby point  .  d is utilized to control the outputting sensitivity of SDW. B. Spectral Euclidean Distance Weight Although SDW is a very vital notion to blur the image, its single use causes a failure at edges, in which the assumption of slow spatial variations does not work. That is, Gaussian low-pass filtering cannot preserve edges. Therefore, the bilateral filter combines both domain filtering and range filtering that respectively consider the closeness and the similarity between a center pixel and its neighboring pixels. Similarity refers to a degree that one pixel can be similar to another. Compared to closeness, similarity, to a large extent, should have both greater capability and higher accuracy to characterize the weight. In bilateral filter, the Euclidean distance between two pixels is used to denote the similarity weight in gray image and color image, and it shows a good blur effect. In hyperspectral image processing, the Euclidean distance is usually involved in many algorithm designation, since it can directly reflect the similarity between two spectrums. For example, the k -nearest neighbor ( k -NN) method takes advantage of the Euclidean distance to robustly capture the local behavior of underlying background distribution [38]. In this work, we still remain the Euclidean distance as one kind of weights, which can be written as:

 1  d ( , x )  2  SEDW (, x)  exp    e    2  e     d e ( , x )    x

(2)

where de (, x) calculates the Euclidean distance between the neighborhood center x and a nearby point  .  e is utilized to control the outputting sensitivity of SEDW. C. Neighboring Background Pixel Selection In hyperspectral anomaly detection field, some background-weighted methods have been proposed to obtain fresher background. These methods compute the weights by using Euclidean distance or spectral angle distance between the neighboring pixel vectors and the background-averaged vector. Two more similar pixel vectors can derive a bigger weight. When the neighboring background pixels are corrupted with anomalies, there is not a much bad impact on the resulting background, since anomalous pixels in the background are imposed on small weights. By contrast, the proposed BFAD outputs weights via the Euclidean distance between a center pixel and its neighboring pixels. If the center pixel is anomalous, then the same anomalies in neighboring pixels will get the largest weights among all. This phenomenon is unfavorable to reconstruct the center pixel. But this unfavorable phenomenon does not occur if the center pixel is background. Therefore, we should think about a method to ensure that the background does not include anomalies when the center pixel is anomalous. Interestingly, the model of dual window is a good choice. As shown in Fig. 1, the dual window is used to separate a local area into two regions, inner window region (IWR) and outer window region (OWR). IWR is called a guarded region. Since anomalous targets are composed of several pixels most of time, IWR makes sense to prevent the anomalies from moving into the background. And in this view, the size of IWR should be equal to or larger than anomalous targets. OWR is employed to model the local background around the target region. Therefore, two types of weights, including SDW and SEDW, respectively utilize the spatial and spectral information in the OWR in this paper.

3

Fig. 1. The model of dual window.

D. Bilateral-Filter-Based Anomaly Detection Two kinds of weights explain the potential similarity between two pixels from the perspectives of spatial distribution and spectral energy. With these weights, the center pixel can be denoted by its surrounding pixels and obtain very similar spectral features to neighboring background. This operation can be regarded as filtering, which is described as N

output (x)  k 1 (x) SDW (i , x)SEDW (i , x)

(3)

i 1

with the normalization: N

k (x)   SDW (i , x)SEDW (i , x)

(4)

i 1

where  i is ith pixel in OWR at the center pixel x , and N is the number of pixels in OWR. If an entire image is inputted into this filter, anomalous points in the original image will be removed, and we can obtain a fresh background image. Fig. 2 illustrates the real original hyperspectral image and the filtered background image. There are three airplanes served as anomalous targets in the original image, but they disappear and are substituted by the neighboring background in the filtered background image after filtering. It is worth noting that our filter still remains a unique favorable feature of the bilateral filter, that is, the edges of the filtered background image are well preserved, which can be seen in Fig. 2(b). If we let I o be the original image and I b be the filtered background image, an anomaly image I a can be obtained by the following function

I a  I o  Ib

(5)

(a) (b) Fig. 2. (a) The 100th band image scene of original image; (b) the 100th band image scene of filtered background image.

It is clear from (5) that both I a and I b are separated from I o . Since each pixel in I b is a background vector, the energies remained in I a will be near zero for the pixels that are background in I o , while for the pixels that are anomalies in I o , their energies in I a will be great. In the anomaly image I a , each pixel vector I ai , which is a row vector, corresponds to the anomaly degree of spectral response of the pixel x i in the original image I o . Hence, the anomaly image I a can be used to detect anomalies according to the anomaly degree of every pixel. In specific, the greater energy in I a means a greater anomaly degree. The anomaly value d i for pixel x i can be derived by

di  I ai I aiT

(6)

4

According to the aforementioned descriptions, the model of bilateral-filter-based anomaly detection for hyperspectral image is presented as follows: 1. Input: a hyperspectral image and parameters: 1) the sizes of dual window; 2) the value of  d ; 3) the value of  e ; 2. Make max-min normalization of the whole hyperspectral image. Suppose X to be a hyperspectral image with M pixels, X  min( X) each of which has L wavebands, the normalization of hyperspectral image is thus defined as: ; max( X)  min( X) 3. For each pixel x i in the normalized hyperspectral image: 1) Construct the new neighboring background falling into the OWR; 2) Compute the two kinds of weights involving SDW and SEDW through (1) and (2); 3) Calculate the weighted average via (3); End for 4. Obtain the anomaly image using (5); 5. Derive anomaly degrees of the anomaly image by (6); Output: detection maps of the entire image. III. HYPERSPECTRAL DATA EXPERIMENTS A.

Description of Hyperspectral Datasets

The first dataset used for experiments is synthesized by two real hyperspectral images. One is the SanDiego Airport (SanDiegoA) hyperspectral dataset, which was obtained by the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS) hyperspectral spectrometer over the area of the SanDiego airport. It contains 224 spectral channels from 0.4 to 2.5m, 126 of which were used for experiments after removing water absorption, low signal-to-noise ratio (SNR) and bad bands. Its spatial resolution is approximately 3.5m. The other one is HyMap data acquired by HyMap hyperspectral remote sensor in 2003 in the Cook City, Montana, USA. This dataset consists of 126 spectral channels from 0.4 to 2.5m as well. It has a size of 280  800 pixels, in which a smaller area with 90  90 pixels (shown in Fig. 3(b)) is segmented as the background image. From SanDiegoA data, four ground materials: gasoline barrel, housed, tree and planes represented by G, H, T, P (shown in Fig. 3(a)) were added into Hymap data in some certain pixels. The abundances of synthetic targets are shown in Fig. 3(c), and the size of targets from left to right is 4  4 , 3  3 , 2  2 and 11 , respectively. Fig. 3(d) presents the ground truth of synthetic data.

(a)

(b)

5

(c)

(d)

Fig. 3. (a) The image scene of SanDiegoA data; (b) the image scene of HyMap data; (c) the image scene of synthetic data; (b) the ground-truth map of synthetic data.

The second dataset is a smaller area with 100  90 pixels segmented from SanDiegoA data (marked by a white rectangle shown in Fig. 3(a)). There are three airplanes as anomalous targets to be detected in the scene, as shown in Fig. 4(a). And its ground truth is given by Fig. 4(b).

(a) (b) Fig. 4. (a) The image scene of SanDiegoA data; (b) the ground-truth map of SanDiegoA data.

The Pavia Center (PaviaC) hyperspectral dataset was collected from the Reflective Optics System Imaging Spectrometer (ROSIS) sensor during a flight over Pavia Center in northern Italy and is available at http://www.ehu.es/ccwintco/index.php/Hyperspectral_Remote_Sensing_ Scenes. The spatial resolution is about 1.3 m per pixel. The image has 115 spectral channels in wavelengths ranging from 430 to 860 nm. After removing the bands that correspond to the water absorption regions, low SNR, and bad bands, 102 available bands of the data are retained in the experiments. In the experiments, a smaller subset with a size of 115 115 pixels, shown in Fig. 5(a), was segmented from the initial larger image. The smaller image constitutes the background, including a bridge, water and shadows, and anomalies representing vehicles on the bridge, which is shown in Fig. 5(b).

(a) (b) Fig. 5. (a) The image scene of PaviaC data; (b) the ground-truth map of PaviaC data.

B.

Experimental Analysis of Parameters

The proposed BFAD has three types of parameters, including the size of the dual window,  d for SDW and  e for SEDW. All of these parameters play significant roles in the resulting detection performance. The receiver operating characteristics (ROC) [39] curve representing detection probability versus false-alarm rates and its area under ROC curve (AUC) [40] are effective techniques to present quantitative performance analysis. A larger AUC represents a better detection performance. The

6

method of global search among four parameters (inner window size, outer window size,  d ,  e ) is used to look for their optimal values, which occurs when the corresponding AUC (inner window size, outer window size,  d ,  e ) is at maximum. Once we obtain their optimal values, three groups of experiments are conducted to investigate their respective effects. TABLE I EXPERIMENTAL PARAMETERS FOR  d Hyperspectral data Synthetic data SanDiegoA data PaviaC data

AND

e

d

e

10 10 10

0.2 1.2 1.0

First, we illustrate the effect of the dual window size on the detection performance of BFAD using three experimental data. The parameter settings of SDW and SEDW are given in Table I. Fig. 6 shows the detection performance about varying dual window sizes. For synthetic data shown in Fig. 6(a), the value of AUC climbs the top when the dual window size is 7/13, but it descends sharply when the inner window size is less than 7, which is expected because anomalous points falling into the OWR can result in a very bad background reconstruction at center pixel. This phenomenon also exists in the detection results of SanDiegoA data and PaviaC data. Their optimum dual window sizes are 13/17 and 9/13, respectively. The experimental results are also the explanation why the dual window is used in this work. Furthermore, it can be concluded that the inner window size is generally larger than that of anomalous targets.

(a)

(b)

7

(c) Fig. 6. Detection performance with respect to the different sizes of dual window for three hyperspectral data. (a) synthetic data; (b) SanDiegoA data; (c) PaviaC data. TABLE II EXPERIMENTAL PARAMETERS FOR DUAL WINDOW SIZE AND  e Hyperspectral data Synthetic data SanDiegoA data PaviaC data

e

Dual window size 7/13 13/17 9/13

0.2 1.2 1.0

Second, we investigate the effect of  d on detection performance of BFAD. Table II presents the parameter selections of dual window size and SEDW. The AUC values with regard to the different  d are given in Table III. For synthetic data, the AUC performance achieves the optimum when  d =5 . For the rest datasets, their largest AUC values are obtained when

 d =10 . All three data basically remain the best detection performance with an increasing value of  d , even if  d =1000 . When  d decreases to 0.1, the corresponding BFAD detects nothing. But on the whole, SDW is insensitive for the output of BFAD. TABLE III EFFECT OF  d

 d =0.1  d =1  d =5  d =10  d =50  d =100  d =1000

ON THE DETECTION PERFORMANCE FOR THREE DATASETS

Area under ROC curve (AUC) Synthetic data SanDiegoA data PaviaC data 0 0 0 0.9996

0.9963

0.9989

0.9997

0.9978

0.9991

0.9997

0.9982

0.9993

0.9997

0.9982

0.9993

0.9997

0.9980

0.9993

0.9997

0.9981

0.9992

Third, we analyze the effect of  e on the detection performance of BFAD. The parameter settings of dual window and SDW are given in Table IV. Fig. 7 presents the ROC curves and the corresponding AUC curves in terms of the changing  e for three groups of hyperspectral data. For synthetic data shown in Fig. 7(a), the detection performance of BFAD upgrades obviously as  e climbs from 0.1 to 0.2. After that, the detection performance gradually declines as  e rises. For SanDiegoA data in Fig. 7(b), the detection performance of BFAD improves as the  e increases to 1.2. And then the detection performance gradually decreases as  e increases. For PaviaC data in Fig.7 (c), BFAD performs best when  e is equal to 1.0, and its detection performance drops gradually when  e increases or decreases. Overall, a small variation of  e will cause a relatively great fluctuation of the detection result. This phenomenon can form a conclusion that BFAD is sensitive to  e .

8

(a)

(b)

(c) Fig. 7. ROC and AUC performance with respect to the different values of  e for three hyperspectral data. (a) synthetic data; (a) SanDiegoA data; (b) PaviaC data; TABLE IV EXPERIMENTAL PARAMETERS FOR DUAL WINDOW SIZE AND  d Hyperspectral data Synthetic data SanDiegoA data PaviaC data

C.

Dual window size 7/13 13/17 9/13

d 10 10 10

Detection Performance

In this section, several experiments are performed to test the detection performance of the proposed BFAD. And experimental parameter settings of BFAD for three datasets are shown in Table V.

9

First, we utilize two real hyperspectral data to illustrate the processing procedures of BFAD step by step. For SanDiegoA data, there are three airplanes in the image shown in Fig. 8(a); but these airplanes disappear from the image and are substituted by the background shown in Fig. 8(b). And in anomaly image scene shown in Fig. 9(c), three airplanes stand out against the background. Such processing effects can be also found for PaviaC data shown in Fig. 9. It should be noted that some anomalies still exist in the background scene of PaviaC as shown in Fig. 9(b). This is because these anomalies are close to each other, causing inevitable background contamination by the anomalies in the OWR.

(a) (b) (c) Fig. 8. Processing procedures of BFAD for SanDiegoA data. (a) 100th band image scene of original image; (b) 100th band image scene of reconstructed background image; (c) 100th band image scene of anomaly image.

(a) (b) (c) Fig. 9. Processing procedures of BFAD for PaviaC data. (a) 100th band image scene of original image; (b) 100th band image scene of reconstructed background image; (c) 100th band image scene of anomaly image. TABLE V EXPERIMENTAL PARAMETERS FOR THREE HYPERSPECTRAL DATA Hyperspectral data Synthetic data SanDiegoA data PaviaC data

Dual window size 7/13 13/17 9/13

d

e

10 10 10

0.2 1.2 1.0

Furthermore, we evaluate the detection performance of our proposed BFAD, comparing with GRXD, LRXD, BACON, RSAD, PAD, LRSMD and MML-RXD. Fig. 10 presents ROC curves of these detectors for three datasets. For synthetic data shown in Fig. 10(a), it is clear that BFAD performs the best among all detectors. When the false alarm rate is a very small value (about 0.05), its probability of detection climbs to 100%. For SanDiegoA data, as shown in Fig. 10(b), the result of BFAD demonstrates a better performance than other anomaly algorithms. And its detection probability is up to 100% with a false alarm rate of around 0.04. Compared with BACON, some improvements can be found in RSAD and PAD. For PaviaC data shown in Fig. 10(c), the ROC curve of BFAD is always above those of other competitors. And when achieving a 100% probability of detection, BFAD has the smallest false alarm rate of approximately 0.004.

10

(a)

(b)

(c) Fig. 10. ROC curves of BFAD and its competitors for three datasets. (a) synthetic data; (b) SanDiegoA data; (c) PaviaC data.

Next, we present two-dimensional (2-D) detection plots of all comparison algorithms with three data sets. Compared with other anomaly detectors, BFAD has higher outputting detection values for anomalous targets among three datasets. For synthetic data as shown in Fig. 11, BFAD and MML-RXD have clearly distinguishable detection outputs, but BFAD can detect more anomalous pixels. For SanDiegoA data and PaviaC data as shown in Fig. 12 and Fig. 13, BFAD also shows this characteristic of distinguishable detection results. It should be noted that BFAD shows a relatively weak capability of background suppression. Since each pixel in the background image is constructed by its neighboring pixels, the difference between it and its corresponding pixel in original image inevitably exists. And this difference can be outputted in the 2-D plot. However, the processing style of BFAD has an advantage in wiping off large blocks of background clutters, which appear in, such as Fig. 11(d), Fig. 12(c) and Fig. 12(f).

(a)

(b)

(c)

(d)

(e) (f) (g) (h) Fig. 11. Two-dimensional plots of the detection results for synthetic data. (a) GRXD; (b) LRXD; (c) BACON; (d) RSAD; (e) PAD; (f) LRSMD; (g) BFAD; (h) MML-RXD.

11

(a)

(b)

(c)

(d)

(e) (f) (g) (h) Fig. 12. Two-dimensional plots of the detection results for SanDiegoA data. (a) GRXD; (b) LRXD; (c) BACON; (d) RSAD; (e) PAD; (f) LRSMD; (g) BFAD; (h) MML-RXD.

For synthetic data shown in Fig. 14(a), the anomaly box is overlapped with the background box for GRXD. The gaps between two boxes for the other seven detectors are very obvious, especially for BFAD that has a slightly bigger gap. The whisker of background box is far from the anomaly box for BFAD, which means our proposed method can effectively suppress the background clutter. For SanDiego data shown in Fig. 14(b), BFAD exhibits a distinct advantage over its competitors in the gap between the anomaly and background boxes. Compared with other detectors, anomaly box of BFAD has the smallest value of the upper whisker that does not go beyond the background box. For PaviaC data shown in Fig. 14(c), all the detectors have obvious gaps between the anomaly and background boxes, and the biggest gap among all comes from BFAD.

(a)

(b)

(c)

(d)

(e) (f) (g) (h) Fig. 13. Two-dimensional plots of the detection results for PaviaC data. (a) GRXD; (b) LRXD; (c) BACON; (d) RSAD; (e) PAD; (f) LRSMD; (g) BFAD; (h) MML-RXD. TABLE VI TIME COMPARISON FOR THREE DATASETS Processing time (seconds) Synthetic data SanDiegoA PaviaC GRXD 0.682 0.708 0.855 LRXD 39.004 42.830 46.180 BACON 169.674 208.236 331.336 RSAD 38.765 58.460 99.668 PAD 12.484 14.368 24.899 LRSMD 2.345 2.981 3.147 BFAD 9.722 13.388 16.581 MML-RXD 37.062 39.926 42.404 TABLE VII COMPUTATIONAL COMPLEXITY OF BFAD SDW

O( N )  O( Nem ) m  R 12

SEDW

O( NL)  O( Nem ) m  R

Also, we report the processing time of the above detectors. All experiments were carried out in MATLAB 2014A on an Intel Core I7-4720 CPU machine with 4 GB of RAM. As shown in Table VI, the processing time of the proposed BFAD is less than that of LRXD, BACON, RSAD, PAD, LRSMD and MML-RXD. And compared with LRXD, BACON, RSAD, PAD and MML-RXD, its computing time decreases by several times or even over ten times. This good detection efficiency can be clearly explained by its computational complexity shown in Table VII. BFAD mainly contains two calculations of SDW and SEDW. Let N denote the number of pixels between outer window and inner window, L denote the number of spectral bands. It is clear from Table VII that the computational complexity of detecting a pixel in BFAD is approximately O( NL)  O( Nem ) where m  R . In contrast, other competing methods usually need to have an inverse operation of L  L matrix when detecting a pixel. Only the computational complexity of this inversion is up to about O( L3 ) , let alone other computations.

(a)

(b)

(c) Fig. 14. ROC curves of BFAD and its competitors for three datasets. (a) synthetic data; (b) SanDiegoA data; (c) PaviaC data.

Finally, we make some discussions about BF, non-local denoising method [41], entropy-based BF [42] and guided BF [35] in terms of hyperspectral anomaly detection. The proposed method is robust to noise. BFAD is aimed to substitute neighboring background pixels for anomalies in a hyperspectral image. After making a subtraction between the original image and the filtered image, pixels that respond to 13

anomalies in the original image will have large values in the resulting image. Based on this operation, our proposed method can easily detect anomalous pixels. After all, noise is very small compared with anomalies, it is hard to damage the whole pixel, especially in high-dimensional hyperspectral image. Therefore, noise involved in hyperspectral image has little effect on the detection. The size of dual window is very important in BFAD. Dual window is used to keep anomalous pixels from participating in the computation of weights, since anomalies can damage the reconstructive background representation when the center pixel is anomalous. Non-local method (NLM) calculates the similarity of two square areas of pixels centered at a neighboring pixel and a center pixel to construct the weight. The non-local means compare the grey level in a single point and the geometrical configuration in a whole neighborhood. However, this method is not suitable for hyperspectral anomaly detection. First, dual window is hard to be utilized in NLM. Second, as mentioned before, the responding weights will be bad when anomalous pixels fall in the neighboring background. In this case, the background cannot be well reconstructed and the detection result is not expected. The key assumption of the guided filter is a local linear model between the guidance I and the filter output q. By this idea, guided filter has the edge-preserving smoothing property. Anomalous targets are usually composed of several pixels, so they cannot be removed through the guided filter like noise. The use of dual window will break the linear model assumption. Therefore, it is difficult to use the guided filter for background reconstruction. Entropy-based BF (EBF) is aimed to become robust to noise by exploiting the information from the denoised estimate and the corresponding method noise, i.e., the difference between a noisy image and its denoised estimate. Moreover, in order to consider the local statistics of images, local entropy is applied to adaptively guide the range parameter selections. As stated before, BFAD is robust to noise in hyperspectral anomaly detection. Therefore, this improved EBF is not very useful for our proposed method. IV. CONCLUSION In this paper, we have proposed a novel hyperspectral anomaly detection algorithm based on the bilateral filter. Different from some existing algorithms that separate the background image from an original image, the proposed algorithm reconstructs the background image according to the original image from a weighting perspective. Inspired by the bilateral filter, we improve it based on the features of hyperspectral image. And the dual window is utilized in two weights to keep anomalous points from falling into the OWR. With these derived weights, the center pixel can be reconstructed as the background via a weighted average of its neighboring pixels. And then anomaly image can be obtained by subtracting the derived background image from the original hyperspectral image. Finally, this anomaly image is used to implement the anomaly detection. Experimental results demonstrate the effectiveness of the proposed BFAD. In the future work, how to adaptively obtain the stable and optimal parameters can be considered. ACKNOWLEDGMENT This work was supported by the National Natural Science Foundation of China (No. 616011351 and No. 61571145), and the Fundamental Research Funds for the Central Universities (No. HEUCF1608). V. REFERENCES D. Landgrebe, “Hyperspectral image data analysis,” IEEE Signal Process. Mag., vol. 19, no. 1, pp. 17–28, Jan. 2002. C.-I. Chang and S.-S. Chiang, “Anomaly detection and classification for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 40, no. 6, pp. 1314–1325, Jun. 2002. [3] L. Zhang, D. Tao, X. Huang, B. Du, and L. Zhang, “Hyperspectral remote sensing image subpixel target detection based on supervised metric learning,” IEEE Trans. Geosci. Remote Sens., vol. 52, no. 8, pp. 4955–4965, Aug. 2014. [4] D. Manolakis and G. Shaw, “Detection algorithms for hyperspectral imaging applications,” IEEE Signal Process. Mag., vol. 19, no. 1, pp. 29–43, Jan. 2002. [5] Y. Yuan, D. Ma and Q. Wang, “Hyperspectral anomaly detection by graph pixel selection,” IEEE Trans. Cybern., vol. 46, no. 12, pp. 3123-3134, Dec. 2015. [6] G. Shaw and D. Manolakis, “Signal processing for hyperspectral image exploitation,” IEEE Signal Process. Mag., vol. 19, no. 1, pp. 12–16, Jan. 2002 [7] S. Matteoli, N. Acito, M. Diani, and G. Corsini, “An automatic approach to adaptive local background estimation and suppression in hyperspectral target detection,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 2, pp. 790–800, Feb. 2011. [8] H. Kwon and N. M. Nasrabadi, “Kernel matched subspace detectors for hyperspectral target detection,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 28, no. 2, pp. 178–194, Feb. 2006. [9] S. Matteoli, T. Veracini, M. Diani and G. Corsini, “Models and Methods for Automated Background Density Estimation in Hyperspectral Anomaly Detection,” IEEE Trans. Geosci. Remote Sens., vol. 51, no. 5, pp. 2837–2852, May. 2013. [10] S. Velasco-Forero, M. Chen, A. Goh and S. Pang, “Comparative Analysis of Covariance Matrix Estimation for Anomaly Detection in Hyperspectral Images,” IEEE J. Sel. Topics Signal Process., vol. 9, no. 6, pp. 1061-1073, Jun. 2015. [1] [2]

14

[11] C.-I. Chang, Hyperspectral Imaging: Techniques for Spectral Detection and Classification. Kluwer Academic/Plenum Publishers: New York, NY, USA, 2003. [12] C.-I. Chang, Hyperspectral Data Processing: Algorithm Design and Analysis. Wiley: Hoboken, NJ, USA, 2013. [13] Q. Guo, B. Zhang, Q. Ran, L. Gao, J. Li, and A. Plaza, “Weighted-RXD and linear filter-based RXD: improving background statistics estimation for anomaly detection in hyperspectral imagery,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 7, no. 6, pp. 2351-2366, Jun. 2014 [14] D. W. J. Stein et al., “Anomaly detection from hyperspectral imagery,” IEEE Signal Process. Mag., vol. 19, no. 1, pp. 58–69, Jan. 2002. [15] C. Zhao, X. Yao, B. Huang, “Real-time anomaly detection based on a fast recursive kernel RX algorithm,” Remote sens., vol. 8, no. 12, Dec. 2016. [16] S. Matteoli, M. Diani, J. Theiler, “An overview of background modeling for detection of targets and anomalies in hyperspectral remotely sensed imagery,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 7, no. 6, pp. 2317-2336, Jun. 2014. [17] N. Nasrabadi, “Hyperspectral Target Detection: An overview of current and future challenges,” IEEE Signal Process. Mag., vol. 31, no. 31, pp. 34-44, Dec. 2013. [18] I. S. Reed and X. Yu, “Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution,” IEEE Trans. Acoust. Speech Signal Process., vol. 38, no. 10, pp. 1760–1770, Oct. 1990. [19] C. Zhao, X. Yao and Y. Yan, “Modified kernel RX algorithm based on background purification and inverse-of-matrix-free calculation,” IEEE Geosci. Remote Sens. Lett., vol. 14, no. 4, pp. 1-5, Feb, 2017. [20] J. Molero, E. Garzon, I. Garcia, A. Plaza, “Analysis and optimizations of global and local versions of the RX algorithm for anomaly detection in hyperspectral data,” IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens., vol. 6, no. 2, pp. 801–814, Apr. 2013. [21] B. Du and L. Zhang, “Random-selection-based anomaly detector for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 5, pp. 1578–1589, May 2011. [22] C.-I. Chang, “Target signature-constrained mixed pixel classification for hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 40, no. 5, pp. 1065–1081, May 2002. [23] N. Billor, A. S. Hadi, and P. F. Velleman, “BACON: blocked adaptive computationally efficient outlier nominators,” Comput. Statist. Data Anal., vol. 34, no. 3, pp. 279–298, Sep. 2000. [24] L. Gao, A. Plaza, B. Zhang, “Probabilistic anomaly detector for remotely sensed hyperspectral data,” J. Appl. Remote Sens., vol. 8, no. 1, Nov. 2014. [25] M. Imani, “RX Anomaly Detector With Rectified Background,” IEEE Geosci. Remote Sens. Lett., vol. 14, no. 8, pp. 1313–1317, Jun. 2017. [26] A. Banerjee, P. Burlina, and C. Diehl, “A support vector method for anomaly detection in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 44, no. 8, pp. 2282–2291, Aug. 2006. [27] Y. Zhang, B. Du, L. Zhang and S. Wang, “A low-rank and sparse matrix decomposition-based mahalanobis distance method for hyperspectral anomaly detection,” IEEE Trans. Geosci. Remote Sens., vol. 53, no. 3, pp. 1376–1389, Mar. 2016. [28] J. Theiler, and G. Grosklos, “Problematic projection to the in-sample subspace for a kernelized anomaly detector,” IEEE Geosci. Remote Sens. Lett., vol. 13, no. 4, pp. 485-489, Apr. 2016. [29] N. Gorelnik, H. Yehudai, and S. R. Rotman, “Anomaly detection in nonstationary backgrounds,” in Proc. 2nd Workshop Hyperspectral Image Signal Process.: Evol. Remote Sens. (WHISPERS 2010), Reykjavik, Iceland, Jun. 14–16, 2010. [30] W. Sun, C. Liu, J. Li, Y. M. Lai and W. Li. “Low-rank and sparse matrix decomposition-based anomaly detection for hyperspectral imagery,” J. Appl. Remote Sens., vol. 8, no. 1, pp. 152-160, May. 2014. [31] M. Golbabaee and P. Vandergheynst, “Hyperspectral image compressed sensing via low rank and joint-sparse matrix recovery,” in IEEE Int. Conf. on Acoustics, Speech and Signal Processing (ICASSP), pp. 2741–2744, 2012. [32] K. Rong et al., “Low-rank and sparse matrix decomposition-based pan sharpening,” in IEEE Int. Geoscience and Remote Sensing Symp. (IGARSS), pp. 2276–2279, 2012. [33] W. Li, Q. Du, “Collaborative representation for hyperspectral anomaly detection,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 3, pp. 1473–1474, Mar. 2015. [34] C. Tomasi and R. Manduchi, “Bilateral Filtering for Gray and Color Images,” International Conference on Computer Vision, pp. 839-846, Jan. 1998. [35] K. He, J. Sun and X. Tang, “Guided image filtering,” IEEE Trans. Pattern Analysis & Machine Intelligence, vol. 35, no. 6, pp. 1397-1409, Jun. 2013. [36] H. Yu, L. Zhao and H. Wang, “Image denoising using trivariate shrinkage filter in the wavelet domain and joint bilateral filter in the spatial domain,” IEEE Trans. Image Process., vol. 18, no. 10, pp. 2364-2369, Jul. 2009. [37] R. Zhao, B. Du, L. Zhang, and L. Zhang, “Beyond background feature extraction: an anomaly detection algorithm inspired by slowly varying signal analysis,” IEEE Trans. Geosci. Remote Sens., vol. 54, no. 3, pp. 1757–1774, Mar. 2016. [38] B. W. Silverman, Density Estimation for Statist. and Data Analysis. London, U.K.: Chapman and Hall, 1986. [39] M. T. Eismann, Hyperspectral Remote Sensing. Bellingham. WA, USA: SPIE Press, 2012, pp. 620–624. [40] S. Khazai, S. Homayouni, A. Safari, and B. Mojaradi, “Anomaly detection in hyperspectral images based on an adaptive support vector method,” IEEE Geosci. Remote Sens. Lett., vol. 8, pp. 646–650, Jul. 2011. [41] A. Buades, C. Bartomeu, and J. M. Morel, “A non-local algorithm for image denoising,” Computer Vision and Pattern Recognition, pp. 60-65, 2005. [42] T. Dai, W. Lu, W. Wang, J. Wang and S. Xia, “Entropy-based bilateral filtering with a new range kernel,” Signal Process., vol. 137, pp: 223-234, Aug. 2017.

15

Highlights

A hyperspectral anomaly detection algorithm based on the bilateral filter is proposed. The background image without anomalous pixels can be obtained via a weighted average. Bilateral-filter-based anomaly detector can effectively suppress the background clutter.

16