literature review

International Journal of Computer Vision 71(2), 127–141, 2007 c© 2006 Springer Science + Business Media, LLC. Manufactured in the United States.

DOI: 10.1007/s11263-006-6689-3

Kernel Spectral Matched Filter for Hyperspectral Imagery

HEESUNG KWON AND NASSER M. NASRABADI Army Research Laboratory, ATTN: AMSRD-ARL-SE-SE, 2800 Powder Mill Rd., Adelphi, MD 20783

[email protected]

[email protected]

Received March 31, 2005; Revised December 12, 2005; Accepted December 14, 2005

First online version published in June, 2006

Abstract. In this paper a kernel-based nonlinear spectral matched filter is introduced for target detection in hyperspectral imagery, which is implemented by using the ideas in kernel-based learning theory. A spectral matched filter is defined in a feature space of high dimensionality, which is implicitly generated by a nonlinear mapping associated with a kernel function. A kernel version of the matched filter is derived by expressing the spectral matched filter in terms of the vector dot products form and replacing each dot product with a kernel function using the so called kernel trick property of the Mercer kernels. The proposed kernel spectral matched filter is equivalent to a nonlinear matched filter in the original input space, which is capable of generating nonlinear decision boundaries. The kernel version of the linear spectral matched filter is implemented and simulation results on hyperspectral imagery show that the kernel spectral matched filter outperforms the conventional linear matched filter.

Keywords: Matched filter· hyperspectral· kernel· nonlinear detection· target detection

1. Introduction

Target detection using linear matched filtering is a well- known approach in detecting objects of interest in hy- perspectral imagery (Manolakis et al., 2000; Robey et al., 1992; Kraut and Scharf, 1999; Kraut et al., 2001; Chang, 2003; Van Veen and Buckley, 1988; Johnson and Dudgeon, 1993; Capon, 1969). Typically, the tar- get spectral signature is obtained from a spectral library or from a set of training data, which is used in con- junction with the correlation (covariance) matrix of the data as the spectral matched filter. However, the linear spectral matched filter detector (MFD) does not ex- ploit the higher order statistical correlations between the spectral bands since it is based only on the sec- ond order statistics. Therefore, its performance is not optimal for non-Gaussian data. Furthermore, the deci- sion boundaries obtained by the conventional MFD is

linear, which is not optimal for non-Gaussian data. The motivation behind designing the nonlinear matched fil- ter is to obtain nonlinear decision boundaries as well as exploiting the higher order statistical correlation be- tween the spectral bands for non-Gaussian data in order to improve the performance of the conventional linear matched filter. Spectral MFDs are based on the assump- tion of a linear model where the spectral signature of the target and the background covariance matrix are as- sumed to be known. A nonlinear spectral matched filter can easily be developed by defining a matched filter like model in a feature space of high dimensionality. How- ever, to implement such a nonlinear matched filter in the feature space may not be computationally tractable due to the high dimensionality of the feature space.

Recently, using the ideas of kernel-based learning theory it has been shown in Müller et al. (2001), Schölkopf et al. (1999), Ruiz and Lopez-de Teruel

128 Kwon and Nasrabadi

(2001), Kwon and Nasrabadi (2004), Baudat and Anouar (2000), Kwon and Nasrabadi (2005) that a number of linear algorithms can easily be extended to nonlinear versions by implementing them in terms of kernel functions, thus avoiding the implementation of the algorithm in the feature space. The input data is in fact implicitly mapped into a kernel feature space where the algorithm is then implemented in terms of certain kernels such as the Mercer kernels (Schölkopf and Smola, 2002).

In this paper, we introduce a kernel-based spectral matched filter in the feature space and drive an expres- sion for its kernel version. We define a matched filter in a kernel feature space, which is equivalent to a non- linear matched filter in the original input space. The matched filter problem is first formulated in a particu- lar kernel feature space, which is implicitly generated by a nonlinear mapping associated with a kernel func- tion. The matched filter expression derived in that fea- ture space is then rewritten in terms of the vector dot products form and by using the so called kernel trick, see (12) in Section 2.2, it is converted in terms of the kernel function. We refer to this process as kernelizing the expression for the nonlinear matched filter and the resulting matched filter is called the kernel-based spec- tral matched filter detector (KMFD). Using the kernel trick idea we avoid implementing the algorithm in the high dimensional feature space. Furthermore, we do not need to know the explicit expression for the nonlinear map that produced the kernel feature space. However, an appropriate kernel function with dot product prop- erty (Mercer kernel) needs to be selected, which has a particular nonlinear mapping associated with it that models the data implicitly in the kernel feature space.

This paper is organized as follows. Section 2 intro- duces the linear matched filter and the idea of kernel trick when using Mercer kernels. In Section 3 the non- linear matched filter is described, which is reformu- lated in terms of the kernel function to obtain the ker- nel matched filter. Performance of the kernel matched filter on hyperspectral imagery is provided in Section 4 and conclusions are given in Section 5.

2. Preliminaries: Linear Matched Filter, Introduction to Kernel Feature Space and Kernel Trick

2.1. Linear Matched Filter

In this section, we introduce the concept of linear spec- tral matched filter. The constrained least squares ap-

proach is used to derive the linear matched filter. Let the input spectral signal x be x = [x(1), x(2), . . . , x(J )]T consisting of J spectral bands. We can model each spectral observation as a linear combination of the tar- get spectral signature and noise

x = as + n, (1)

where a is an attenuation constant (target abundance measure). When a = 0 no target is present and when a > 0 target is present, vector s = [s(1), s(2), . . . , s(J )]T contains the spectral signature of the target and vector n contains the additive background clutter noise.

We can design a linear matched filter w = [w(1), w(2), . . . , w(J )]T such that the desired target sig- nal s is passed through while the average filter out- put energy is minimized. Define X to be a J × N matrix of the N reference pixels obtained from the test input image. Let each observation spectral pixel to be represented as a column in the sample matrix X

X = [x1 x2 · · · xN ]. (2)

The output of the filter for the input xi is given by

y(xi ) = wT xi = xTi w. (3)

The average output power of the filter for the reference data X is given by

1

N

N∑ i=1

y(xi )2 = wT (

1

N

N∑ i=1

xi xTi

) w = wT R̂w, (4)

where the R̂ is the estimated correlation matrix of the reference data. This constrained filter design is equiv- alent to a constrained least squares minimization prob- lem, as was shown in Scharf (1991), Van Veen and Buckley (1988), Harsanyi (1993), Chang (2003), which is given by

min w

{wT R̂w} subject to sT w = 1 (5)

where minimization of minw{wT R̂w} ensures the back- ground clutter noise is suppressed by the filter w and the constrain condition sT w = 1 makes sure that the filter gives an output of unity when a target is detected.

Kernel Spectral Matched Filter for Hyperspectral Imagery 129

The solution to this quadratic minimization problem is given by

w = R̂ −1s

sT R̂−1s (6)

where correlation matrix R̂ is usually estimated from the input image. The expression (6) is referred to as the minimum variance distortionless response (MVDR) beamformer in the array processing litera- ture (Van Veen and Buckley, 1988; Johnson and Dud- geon, 1993). More recently, the same expression was derived in the hyperspectral target detection literature and was called the constrained energy minimization (CEM) filter or the spectral correlation-based matched filter (Harsanyi, 1993; Chang, 2003). The output of the linear correlation-based matched filter for a test input r, given the estimated correlation matrix, is given by

y(r) = wT r = s T R̂−1r

sT R̂−1s . (7)

If the mean of the observation data is removed (cen- tered) a similar expression is obtained for the centered data which is given by

yr = wT r = s T Ĉ−1r

sT Ĉ−1s (8)

where Ĉ represents the estimated covariance matrix for the centered reference image data. In the above derivation of the correlation-based matched filter no assumption was made about the distribution of the ad- ditive noise n in the linear model (1). However, if n is assumed to be a Gaussian random noise distributed as N (0, C), it has been shown in Robey et al. (1992) and Kraut and Scharf (1999) that using the Generalized Likelihood Ratio Test (GLRT) a similar expression to (8), as in MVDR or CEM, can be obtained for the es- timated abundance measure â given by

â = s T C−1r

sT C−1s . (9)

It should be noted that C is now the expected covari- ance matrix of the background noise only and does not include any target data. When the estimated covari- ance matrix Ĉ of the background data is used in (9) this filter is referred to as the adaptive matched filter (Robey et al., 1992) in the signal processing literature or Capon method (Capon, 1969) in the array processing

literature. In Robey et al. (1992) it was shown that the CFAR behavior of this filter is given by

α(r) = |s T Ĉ−1r|2 sT Ĉ−1s

(10)

which is proportional to the estimated squared magni- tude of the output matched filter referred in Kraut et al. (2001) as the signal-to-noise ratio (SNR).

In this paper, we only present the experimental re- sults for the linear correlation-based matched filter given by the expression (7). Similar results are obtained by using (8) for the centered data.

2.2. Kernel Feature Space and Kernel Trick

In this subsection an introduction to kernel feature map and kernel learning is provided, which is used in the next section to convert a nonlinear version of the matched filter into its corresponding kernel format. Suppose the input hyperspectral data is represented by the data space (X ⊆ RJ ) and F be a nonlinear fea- ture space associated with X by a nonlinear mapping function φ

φ : X → F, x �→ φ(x), (11)

where x is an input vector in X , which is mapped into a potentially much higher dimensional feature space.

Any linear algorithm can be remodeled in this high dimensional feature space by replacing the original in- put data x with the mapped data φ(x). Implementing any linear algorithm (e.g., matched filter) in the feature space is equivalent to performing a nonlinear version of that algorithm (i.e., nonlinear matched filter) in the original data space. Due to the high dimensionality of the feature spaceF it is computationally not feasible to implement the algorithm in the feature space. However, in the kernel-based learning algorithms the task is first formulated in terms of dot products in the feature space and then the kernel trick (12) is used to implement the dot products in terms of kernel functions (Schölkopf and Smola, 2002). The kernel representation for the dot products in F (known as the kernel trick) is expressed as

k(xi , x j ) = 〈φ(xi ), φ(x j )〉 = φ(xi ) · φ(x j ), (12)

130 Kwon and Nasrabadi

where k is a positive definite kernel, such as a Mercer kernel (Schölkopf and Smola, 2002).

Using this kernel trick allows us to implicitly com- pute the dot products in F without mapping the input vectors into F ; therefore, in the kernel-based learn- ing methods, the mapping φ does not need to be iden- tified. However, an appropriate kernel has to be de- fined, which has a nonlinear mapping associated with it. Equation (12) shows that the dot products in F can be avoided and replaced with a kernel, which can be easily calculated without identifying the nonlinear map φ. Experimental results for three different kernels are reported in this paper. The three kernels are

(i) the Gaussian Radial Bases Function kernel (RBF kernel): k(x, y) = exp(−‖x−y‖2c ),

(ii) spectral angle-based kernel: k(x, y) = x·y‖x‖‖y‖ , and

(iii) polynomial kernel: k(x, y) = ((x · y) + θ )d

where c represents the width of the Gaussian RBF ker- nel, d is a positive integer, and θ ≥ 0 is a constant. See Schölkopf and Smola (2002) for more detailed infor- mation about the properties of kernels and kernel-based learning theory.

3. Nonlinear Matched Filter and Kernel Matched Filter

In this section, we show how to formulate a kernel version of the linear matched filter. This is achieved by modeling the data in a kernel feature space where the corresponding linear matched filter in the feature space is equivalent to a nonlinear matched filter in the input space. We then show how to implement the nonlinear matched filter in terms of the kernel function.

3.1. Introduction to Nonlinear Matched Filter

Consider a linear model of the input data in a kernel feature space given by

φ(x) = aφφ(s) + nφ, (13)

where φ is a nonlinear mapping associated with a ker- nel function, aφ is an attenuation constant (target abun- dance measure), the high dimensional vector φ(s) con- tains the spectral signature of the target in the feature space, and vector nφ contains the additive noise in the

feature space. The above linear model (13) in the fea- ture space is not the same as the nonlinearly mapped version of the additive model given in (1). However, this linear model in the feature space is equivalent to a specific nonlinear model in the input space. Therefore, defining a matched filter using the linear model (13) is the same as developing a nonlinear matched filter for a specific nonlinear model in the input space.

Using the constrained least squares approach that was explained in the previous section it can easily be shown that the equivalent correlation-based matched filter wφ in the feature space is given by

wφ = R̂−1φ φ(s)

φ(s)T R̂−1φ φ(s) , (14)

where R̂φ is the estimated correlation of pixels in the feature space. The estimated correlation matrix is given by

R̂φ = 1 N

XφXφT (15)

where Xφ = [φ(x1) φ(x2) . . . φ(xN )] is a matrix whose columns are the mapped input reference data in the feature space. The matched filter in the feature space (14) is equivalent to a nonlinear matched filter in the input space and its output for the input φ(r) is given by

y(φ(r)) = wTφ φ(r) = φ(s)T R̂−1φ φ(r)

φ(s)T R̂−1φ φ(s) . (16)

Due to the high dimensionality of the feature space the expression (16) is not tractable. Therefore, we can- not directly implement it in the feature space. We need to first convert this expression in terms of the dot prod- ucts of the input vectors in the feature space. The kernel trick is then used to convert the dot products in the fea- ture space in terms of the kernel functions.

3.2. Kernel Matched Filter

In this subsection, we show how to kernelize the matched filter in the feature space. The estimated back- ground correlation matrix can be represented by its eigenvalue decomposition or so called spectral decom- position (Strang, 1986) given by

R̂φ = Vφ�VφT , (17)

Kernel Spectral Matched Filter for Hyperspectral Imagery 131

Figure 1. Two-dimensional toy data sets: (a) A Gaussian mixture and (b) nonlinearly mapped data.

Figure 2. Contour and surface plots of MSD and KMSD: (a) MFD on the data shown in Fig. 1(a), (b) KMFD on the data shown in Fig. 1(a), (c) MFD on the data shown in Figs. 1(b) and (d) KMFD on the data shown in Fig. 1(b).

where � is a diagonal matrix consisting of the nonzero eigenvalues of R̂φ and Vφ is a matrix whose columns are the eigenvectors of R̂φ in the feature space. The eigenvector matrix is represented by

Vφ = [v1φ v2φ · · · vNφ ], (18)

where N is the maximum number of eigenvectors with nonzero eigenvalue.

The sample correlation (covariance) matrix in the feature space is rank deficient because the number of samples are usually much less than the dimensional- ity of the feature space. Therefore, the inverse of R̂φ

cannot be obtained and we have to resort to the pseudo- inverse of the sample correlation matrix, which is the minimum length least squares solution for the inverse sample correlation matrix (see p. 450 in Strang, 1986). The pseudo-inverse of the estimated background cor- relation matrix can be written in terms of its eigen- value decomposition or singular value decomposition (Strang, 1986) as

R̂#φ = Vφ�−1VφT . (19)

Furthermore, the diagonal matrix �−1 can be replaced with a truncated version of �−1 by only including

132 Kwon and Nasrabadi

Figure 3. Sample band images (48th) from the HYDICE and mine images. (a) Desert Radiance II image, (b) Forest Radiance I image and (c) mine image.

the eigenvalues that are above a small threshold in order to obtain what is called the effective rank of the matrix (see p. 445 in Strang, 1986). Truncation of the diagonal matrix �−1 will numerically provide a more stable pseudo-inverse since, due to round-off errors, it is not easy to identify the true non-zero eigenvalues.

Each eigenvector v jφ in the feature space, as shown in Schölkopf et al. (1999), can be expressed as a linear combination of the input reference vectors in the feature space given by

v jφ = N∑

i=1 λ j

−1/2β ji φ(xi ) = Xφβ jλ j −1/2, (20)

where the expansion coefficient vectors β j = (β j1 , β

j 2 , . . . , β

j N )

T , for j = 1, . . . , N1, N1 ≤ N , are the eigenvectors with nonzero eigenvalues of the ker- nel (Gram) matrix K(X, X) as shown in Appendix I,

K(X, X) = (K)i j is an N × N matrix whose entries are the dot products k(xi , x j ) =< φ(xi ), φ(x j ) > for xi , x j ∈ X; and λ j , j = 1, . . . , N1, N1 ≤ N are the corresponding nonzero eigenvalues associated with the eigenvectors. For all the eigenvectors Vφ in the feature space we have

Vφ = XφB�−1/2, (21)

where B = [β1 β2 . . . βN1 ]. Substituting (21) into (19) yields

R̂#φ = XφB�−2BT XTφ . (22)

Inserting Eq. (22) into (16) it can be rewritten as

y(φ(r)) = φ(s) T XφB�−2BT XTφ φ(r)

φ(s)T XφB�−2BT XTφ φ(s) . (23)

Kernel Spectral Matched Filter for Hyperspectral Imagery 133

Figure 4. Examples of spectral curves from various regions in the Forest Radiance I image: (a) tree, (b) grass, (c) shadow, and (d) target regions. In general, the spectral curves show a wide range of spectral variability. Especially, the target spectral curves are remarkably different in a local region.

The dot product term φ(s)T Xφ in the feature space can be represented in terms of the kernel function, which is referred to as its empirical kernel map (Schölkopf and Smola, 2002)

φ(s)T Xφ = (k(x1, s), k(x2, s), . . . , k(xN , s)) = kT (X, s) = kTs . (24)

Similarly ,

φ(r)T Xφ = (k(x1, r), k(x2, r), . . . , k(xN , r)) = kT (X, r) = kTr . (25)

Also using the properties of the Kernel Principal Com- ponent Analysis (PCA) as shown in Appendix I, we have the relationship

K−2 = 1 N 2

B�−2BT . (26)

Substituting (24), (25), and (26) into (23) the kernelized version of the matched filter is given by

y(kr) = k(X, s) T K−2k(X, r)

k(X, s)T K−2k(X, s) = k

T s K

−2kr kTs K−2ks

, (27)

which can now be implemented with no knowledge of the mapping function φ. The only requirement is a good choice for the kernel function k.

In the expression (27) the inverse K−2 may not be numerically stable if the background spectral samples are not independent. Therefore, the pseudo-inverse of K is used, which is based on eigenvalue decomposition in (42) where eigenvectors with non-zero eigenvalues are used. In the experimental section, expression (42) is used to obtain the pseudo-inverse of K where only the eigenvectors with eigenvalues above a small threshold are kept. A similar procedure has been used in Ruiz and Lopez-de Teruel (2001) to obtain a stable pseudo- inverse of the Gram matrix by discarding the lowest eigenvalues that are below 10−5. The number of eigen- vectors that are kept will determine the effective rank of the matrix K.

In the derivation of the correlation-based kernel matched filter we assumed that the data was not cen- tered in the feature space. To obtain covariance-based kernel matched filter, expression (8), we need to cen- ter the data by removing the sample mean in the feature space. However, removing the sample mean is not computationally tractable in the feature space due to the high dimensionality of F . Therefore, the

134 Kwon and Nasrabadi

Figure 5. Examples of spectral curves from various regions in the Desert Radiance II image: (a) vegetation, (b) dirt road, (c) soil, and (d) target regions. The target spectral curves are also characterized by a wide range of spectral variability.

kernel matrix K needs to be properly centered, as shown in Schölkopf and Smola (2002). The effect of center- ing on the kernel matched filter can be incorporated by replacing the uncentered Xφ with the centered Xφ −μφ (where μφ = 1N

∑N i=1 φ(xi ) is the mean of the refer-

ence input data) in the estimation of the centered cor- relation (covariance) matrix expression (15) as well as (24) and (25) for the empirical kernel mappings of the target and input data, respectively. The resulting cen- tered K̂ is shown in Schölkopf and Smola (2002) to be given by

K̂ = (K − 1N K − K1N + 1N K1N ), (28)

where the elements of the N ×N matrix (1N )i j = 1/N . The properly centered kernel matched filter output for (27) corresponding to the kernel version of (8) is now given by

y(k̂r) = k̂ T s K̂

−2k̂r k̂Ts K̂−2k̂s

, (29)

where k̂Ts = kTs − 1N ∑N

i=1 k(xi , s)�1 and k̂Tr = kTr − 1 N

∑N i=1 k(xi , r)�1, which are obtained by replacing Xφ

with Xφ − μφ in (24), and (25), respectively. The col-

umn vector �1 denotes an N-dimensional vector with all its components equal to 1.

Furthermore, the CFAR behavior of the kernel matched filter (29) is given by

α(k̂r) = |k̂ T s K̂

−2k̂r|2 k̂Ts K̂−2k̂s

. (30)

4. Simulation Results

In this section, we implemented both the proposed KMFD described by (27) and the conventional MFD described by (10) on simulated data as well as real hyperspectral imagery. We implemented KMFD with three kernel functions, each kernel function being as- sociated with a different feature space. The three dif- ferent kernels used were (i) the Gaussian RBF kernel, exp(−‖x−y‖

2

30 ), (ii) spectral angle-based kernel, x·y

‖x‖‖y‖ , and (iii) 5th order polynomial kernel, ((x · y) + 1)5.

The simulated data consists of two illustrative two- dimensional toy data sets; a Gaussian mixture, as shown in Fig. 1(a), and a nonlinearly mapped data shown in Fig. 1(b). In Fig. 1 data points for the desired target were represented by the star-shaped symbol and

Kernel Spectral Matched Filter for Hyperspectral Imagery 135

Figure 6. Detection results for the Desert Radiance II image using the kernel matched filter detectors (KMFD) and the matched filter detector (MFD): (a) KMFD with the Gaussian RBF kernel, (b) KMFD with the polynomial kernel, (c) KMFD with the spectral angle-based kernel, and (d) MFD in the original input domain.

the background points were represented by the circles. In Fig. 1(b) the two-dimensional data points x = (x, y) for each class was obtained by nonlinearly mapping the original Gaussian mixture data points x0 = (x0, y0) in Fig. 1(a). All the data points in Fig. 1(a) were nonlin- early mapped by x = (x, y) = (x0, x20 + y0), therefore, the second component of each data point is nonlinearly related to its first component.

Figure 2 shows contour and surface plots of the con- ventional MFD and KMFD on the toy data sets as

shown in Fig. 1. For both data sets, the contours gen- erated by KMFD are highly nonlinear and naturally follow the dispersion of the data, thus successfully sep- arating the two classes, as opposed to the linear con- tours obtained by MSD. For both the Gaussian mixture and nonlinearly mapped data, KMFD clearly provided significantly improved discrimination over the conven- tional MFD.

The real hyperspectral images are from a HY- DICE (HYperspectral Digital Imagery Collection

136 Kwon and Nasrabadi

Figure 7. Detection results for the Forest Radiance I image using the kernel matched filter detectors (KMFDs) and the matched filter detector (MFD): (a) KMFD with the Gaussian RBF kernel, (b) KMFD with the polynomial kernel, (c) KMFD with the spectral angle-based kernel and (d) MFD in the original input domain.

Experiment) sensor and a mine detection sensor. The HYDICE imaging sensor generates 210 bands across the whole spectral range (0.4–2.5 μm). But we only use 150 bands by discarding water absorption and low SNR bands; the spectral bands used are the 23rd– 101st, 109th–136th, and 152nd–194th. The hyperspec- tral mine image consists of 70 bands over the spectral range of 8–11.5 μm which includes the long-wave in- frared (LWIR) band. Two HYDICE images from the Desert Radiance II data collection and the Forest Ra-

diance I data collection, and the mine image were used to test both the kernel-based and conventional matched filter detectors. The Desert Radiance II (DR-II) image contains 6 military targets located in the dirt road; the Forest Radiance I (FR-I) image includes total of 14 military targets along the tree line; and the hyperspec- tral mine image contains a total of 33 surface mines, as shown in the sample band images in Fig. 3. While spectral curves from various regions of different mate- rials (or terrain types) show somewhat different spectral

Kernel Spectral Matched Filter for Hyperspectral Imagery 137

Figure 8. Detection results for the mine image using the kernel matched filter detectors (KMFDs) and the matched filter detector (MFD): (a) KMFD with the Gaussian RBF kernel, (b) KMFD with the polynomial kernel kernel, (c) KMFD with the spectral angle-based kernel kernel and (d) MFD in the original input domain.

characteristics, even in a local region of the same ma- terial types they show a wide range of spectral variabil- ity (especially for targets and vegetation), thus making target detection a challenging task, as shown in Figs. 4 and 5.

All the pixel vectors in a test image are first normal- ized by a constant, which is a maximum value obtained from all the spectral components of the spectral vec- tors in the corresponding test image, so that the entries of the normalized pixel vectors fit into the interval of spectral values between zero and one. The rescaling of

pixel vectors was mainly performed to effectively uti- lize the dynamic range of Gaussian RBF kernel. The rescaling does not affect the performance of KMFDs or MFD when the other two kernel functions are used.

4.1. Algorithm Implementation

Background statistics—the kernel matrix K̂ and the correlation matrix R̂—are globally estimated for the kernel-based and conventional implementation of the

138 Kwon and Nasrabadi

Figure 9. ROC curves obtained from the detection results for the Desert Radiance II image shown in Fig. 6.

Figure 10. ROC curves obtained from the detection results for the Forest Radiance I image shown in Fig. 7.

matched filter detectors, respectively. Global estima- tion must be performed prior to detection and normally needs a large amount of data samples to successfully represent all the background types present in a given data set. In this paper, we use a large number of spectral

vectors obtained from a given test image that best rep- resent the spectral characteristics of the background. A well-known data clustering algorithm, k-means (Jain et al., 1999), is used on the spectral vectors in order to generate a significantly less number of spectral vectors

Kernel Spectral Matched Filter for Hyperspectral Imagery 139

Figure 11. ROC curves obtained from the detection results for the mine image shown in Fig. 8.

(centroids) from which appropriate background statis- tics are estimated. The number of the representative spectral vectors obtained from the k-means procedure was set to 300 in the experiment. The correlation and kernel matrices were obtained from these 300 spectral vectors. In all the experiments the same target spectral signature s is obtained by averaging the target samples collected from one of the targets in the test image. For both the DR-II and FR-I images it is the left most target.

Since a relatively small set of the background sam- ples was used to represent the background statistics of the given test images, the kernel matrix K̂ and the cor- relation matrix R̂ need to be regularized in order to obtain a numerically stable pseudo-inverse. In this pa- per, in order to regularize both K̂ and R̂ the eigenvalues decomposition of each of the two matrices is first per- formed. The eigenvectors with large eigenvalues con- vey more relevant information to the prior knowledge of the background than the ones with very small eigen- values that could represent sensor noise. The regular- ized versions of K̂ and R̂ are obtained by discarding all the eigenvalues below a small threshold, 10−5. In our experimental results we did not have to regularize the correlation matrix since the number of the back- ground data was sufficient to obtain the correct inverse. However, in the case of inverse Gram matrix the small eigenvalues were always discarded.

4.2. Performance Comparison

The receiver operating characteristics (ROC) curves representing detection probability Pd versus false alarm rates N f were generated to provide quantitative performance comparison as well as qualitative perfor- mance comparison. For ROC curves generation, based on the ground truth information for the HYDICE im- ages, we obtain the coordinates of all the rectangular target regions. Each target was considered to be de- tected if at least one pixel within the corresponding target region was detected at a given false alarm rate. Pd and N f are defined as

Pd := NhitNt and N f := Nmiss Ntot

, (31)

where Nhit represents the number of targets detected given a certain threshold; Nt represents the total num- ber of targets in the hyperspectral image; Nmiss rep- resents the number of background pixels detected as targets; and Ntot represents the total number of pixels in the hyperspectral image. Pd becomes one only when all the targets are detected.

Figures 6–8 show the detection results for the DR- II, FR-I and mine images using KMFD with the three different kernels and correlation-based MFD, respec- tively. The corresponding ROC curves for the detection

140 Kwon and Nasrabadi

results are shown in Fig. 9–11. For the DR-II image, as shown in Fig. 9, KMFD with any choice of the three kernels clearly outperformed the conventional MFD at almost every false alarm rate. For the FR-I image the background structure is much more complex than that of DR-I. It includes the tree area where most irregular illumination effects occur, the long shadowy transition and the region filled mostly with grass. KMFD using any choice of the kernels still showed improved de- tection results over the conventional MFD for both the FR-I and the mine images, as shown in Figs. 10 and 11, respectively.

5. Conclusions

We have extended the conventional matched filter de- tector to a nonlinear version by implicitly mapping the input data into a much higher dimensional feature space in order to make use of high-order nonlinear correla- tions between the spectral bands of a hyperspectral im- age. The expression of the nonlinear matched filter de- tector in the feature space, which is basically intractable due to high (potentially infinite) dimensionality, is con- verted in terms of kernels using the kernel eigenvector representation as well as the kernel trick to derive a tractable algorithmic expression.

KMFD, the kernel counterpart of MFD, was im- plemented with several different kernels, each with different characteristics. In general, KMFD with all the kernels showed a superior detection performance when compared to the conventional MFD for the HY- DICE and mine images tested in this paper. The detec- tion results show that the kernel-based nonlinear detec- tion method is quite suitable for identifying underlying structures of complex data such as hyperspectral data, thus they are more powerful in discriminating targets of interest.

Appendix I: Kernel PCA

In this Appendix, we present derivation of Kernel PCA and its properties providing the relationship between the correlation (covariance) matrix and the correspond- ing Gram (centered) matrix. Our goal is to prove (26). To drive the Kernel PCA consider the estimated back- ground clutter correlation matrix in the feature space and assume that the input data is not centered . The estimated correlation matrix in the feature space is

given by

R̂φ = 1 N

XφXTφ . (32)

The PCA eigenvectors are computed by solving the eigenvalue problem

λvφ = R̂φvφ

= 1 N

N∑ i=1

φ(xi )φ(xi )T vφ

= 1 N

N∑ i=1

〈φ(xi ), vφ〉φ(xi ), (33)

where vφ is an eigenvector in F with a corresponding nonzero eigenvalue λ. Equation (33) indicates that each eigenvector vφ with corresponding λ �= 0 are spanned by φ(x1), . . . , φ(xN )—i.e.

vφ = N∑

i=1 λ−1/2βiφ(xi ) = λ−1/2Xφβ, (34)

where Xφ = [ φ(x1) φ(x2) . . . φ(xN ) ] and β = (β1, β2, . . . , βN )T . Substituting (34) into (33) and mul- tiplying with φ(xn)T , n = 1, . . . , N , yields

λ

N∑ i=1

βi < φ(xn), φ(xi ) >

= 1 N

N∑ i=1

βiφ(xn)φ(xi )φ(xi )T N∑

i=1 φ(xi )

= 1 N

N∑ i=1

βi 〈φ(xn), N∑

j=1 φ(x j )〈φ(x j ), φ(xi )〉〉,

for all n = 1, . . . , N . (35)

We denote by K = K(X, X) = (K)i j the N × N kernel matrix whose entries are the dot products 〈φ(xi ), φ(x j )〉. Equation (35) can now be rewritten as

Nλβ = Kβ, (36)

where β turn out to be the eigenvectors with nonzero eigenvalues of the kernel matrix K. Therefore, the Gram matrix can be written in terms of its eigenvector decomposition as

K = B�BT , (37)

Kernel Spectral Matched Filter for Hyperspectral Imagery 141

where B = [β1 β2 . . . βN ] are the eigenvectors of the kernel matrix and � is a diagonal matrix with diagonal values equal to the eigenvalues of the kernel matrix K. Similarly, from the definition of PCA in the feature space (33) the estimated background correlation matrix is decomposed as

R̂φ = Vφ�VφT , (38)

where Vφ = [v1φ v2φ . . . vNφ ] and � is a diagonal ma- trix with its diagonal elements being the eigenvalues of Ĉφ . From (33) and (36) the eigenvalues of the correla- tion matrix � in the feature space and the eigenvalues of the kernel matrix � are related by

� = 1 N

�. (39)

Substituting (39) into (37) we obtain the relationship

K = NB�BT , (40)

where N is a constant representing the total num- ber of background clutter samples, which can be ignored.

The sample correlation matrix in the feature space is rank deficient, therefore, its inverse cannot be obtained but its pseudo-inverse can be written as (Strang, 1986)

R̂#φ = Vφ�−1VφT = XφB�−2BT XφT . (41)

The maximum number of eigenvectors in the pseudo- inverse is equal to the number of non-zero eigenvalues (or the number of independent data samples), which cannot be exactly determined due to round-off er- ror in the calculations. Therefore, the effective rank (Strang, 1986) is determined by only including the eigenvalues that are above a small threshold. Simi- larly, the inverse Gram matrix K−1 can also be written as

K−1 = 1 N B�−1BT . (42)

If the data samples are not independent then the pseudo- inverse of the Gram matrix has to be used, which is the same as (42) except only the eigenvectors with eigen- values above a small threshold are included in order to obtain a numerically stable inverse. Using (42) it is obvious that the squared inverse Gram matrix can also

be written as

K−2 = 1 N 2

B�−2BT . (43)

Acknowledgment

The authors would like to thank the anonymous review- ers for their insightful comments and suggestions that helped us to improve the quality of our paper.

References

Baudat, G. and Anouar, F. 2000. Generalized discriminant analysis using a kernel approach. Neural Computat., 12:2385–2404.

Capon, J. 1969. High-resolution frequency-wavenumber spectrum analysis. Proc. of the IEEE, 57:1408–1418.

Chang, C.-I. 2003. Hyperspectral Imaging: Techniques for Detection and Classification. Kluwer Academic/Plenum Publishers.

Harsanyi, J.C. 1993. Detection and Classification of Subpixel Spec- tral Signatures in Hyperspectral Image Sequences. Ph.D. disser- tation, Dept. Elect. Eng., Univ. of Maryland, Baltimore County.

Jain, A.K., Murty, M.N., and Flynn, P.J. 1999. Data clustering: A review. ACM Computing Surveys, 31(3):264–323.

Johnson, D.H. and Dudgeon, D.E. 1993. Array Signal Processing. Prentice Hall.

Kraut, S. and Scharf, L.L. 1999. The CFAR adaptive subspace detec- tor is a scale invariant-invariant GLRT. IEEE Trans. Signal Pro- cess., 47(9):2538–2541.

Kraut, S., Scharf, L.L., and McWhorter, T. 2001. Adaptive subspace detectors. IEEE Trans. Signal Process., 49(1):208–216.

Kwon, H. and Nasrabadi, N.M. 2004. Kernel-based subpixel target detection in hyperspectral images. In Proc. of IEEE Joint Confer- ence on Neural Networks, Budapest, Hungary, pp. 717–722.

Kwon, H. and Nasrabadi, N.M. 2005. Kernel RX-algorithm: A non- linear anomaly detector for hyperspectral imagery. IEEE Trans. Geosci. Remote Sensing, 43(2):388–397.

Manolakis, D., Shaw, G., and Keshava, N. 2000. Comparative analy- sis of hyperspectral adaptive matched filter detector. In Proc. SPIE, vol. 4049, pp. 2–17.

Müller, K.-R., Mika, S., Rätsch, G., Tsuda, K., and Schölkopf, B. 2001. An introduction to kernel-based learning algorithms. IEEE Trans. Neural Networks., 2:181–202.

Robey, F.C., Fuhrmann, D.R., and Kelly, E.J. 1992. A CFAR adaptive matched filter detector. IEEE Trans. on Aerospace and Elect. Syst., 28(1):208–216.

Ruiz, A. and Lopez-de Teruel, E. 2001. Nonlinear kernel-based sta- tistical patten analysis. IEEE Trans. Neural Networks., 12:16–32.

Scharf, L.L. 1991. Statistical Signal Processing. Addison-Wesley. Schölkopf, B. and Smola, A.J. 2002. Learning with Kernels. The

MIT Press. Schölkopf, B., Smola, A.J., and Müller, K.-R. 1999. Kernel principal

component analysis. Neural Computat., (10):1299–1319 . Strang, G. 1986. Linear Algebra and Its Applications. Harcourt Brace

& Company. Van Veen, B.D. and Buckley, K.M. 1988. Beamforming: A versatile

approach to spatial filtering. IEEE ASSP Magazine, pp. 4–24.