An Algorithm for Fully Constrained Abundance ...

1 downloads 0 Views 369KB Size Report
Samuel Rosario-Torres and Miguel Vélez-Reyes. Laboratory for Applied Remote Sensing and Image Processing,. University of Puerto Rico at Mayagüez, P.O. ...
An Algorithm for Fully Constrained Abundance Estimation in Hyperspectral Unmixing Samuel Rosario-Torres and Miguel Vélez-Reyes Laboratory for Applied Remote Sensing and Image Processing, University of Puerto Rico at Mayagüez, P.O. Box 9048, Mayagüez, PR 00681-9048 e-mail: [email protected], [email protected] ABSTRACT This paper presents an algorithm for abundance estimation in hyperspectral imagery. The fully constrained abundance estimation problem where the positivity and the sum to less than or equal to one (or sum equal to one) constraints are enforced is solved by reformulating the problem as a least distance (LSD) least squares (LS) problem. The advantage of reformulating the problem as a least distance problem is that the resulting LSD problem can be solved using a duality theory using a nonnegative LS problem (NNLS). The NNLS problem can then be solved using Hanson and Lawson algorithm or one of several multiplicative iterative algorithms presented in the literature. The paper presents the derivation of the algorithm and a comparison to other approaches described in the literature. Application to HYPERION image taken over La Parguera, Puerto Rico is presented. Keywords: Unmixing, Hyperspectral Imagery

1

INTRODUCTION

Hyperspectral sensors collect hundreds of narrow and contiguously spaced spectral bands of data organized in the socalled hyperspectral cube. Hyperspectral imagery provides fully registered spatial and high-resolution spectral information that is invaluable in discriminating between objects since they have unique spectral signatures that are captured by the data. The spatial resolution of most HSI flown nowadays is larger than the size of the objects being observed; in this case the high spectral sensor resolution can be used to detect and classify subpixel objects by their contribution to the measured spectral signal. The problem of interest is to decompose the measured reflectance (or radiance) into its basic elements. This is the so-called unmixing problem [1] in HSI. Spectral unmixing is the procedure by which the measured spectrum of a pixel is decomposed into a collection of constituent spectra, or endmembers, and a set of corresponding fractions or abundances. This paper presents an algorithm that solves the fully constrained abundance estimation problem. Experimental results using a Hyperion image are presented to demonstrate the algorithm performance.

2

MIXING MODEL

Analytical models for the mixing of materials provide the basis for the development of techniques to recover estimates of the constituents and their proportions from mixed pixels. In the typical mixing process, the surface is portrayed as a checkerboard mixture. The incident light interacts with the surface and the received light conveys the characteristics of the media in a proportion equal to the area covered by the endmembers and the reflectivity of the media. The typical model for such interaction is given by [1] b=

n

∑x a

i i

+ w = Ax + w

(1)

i =1

where b ∈ ℜ m+ is the pixel of interest, ai is the spectral signature of the i-th endmember, xi is the corresponding fractional abundance, w is the measurement noise, m is the number of spectral channel, and n is the number of n is the matrix of endmembers and x ∈ ℜ n+ is the vector of spectral abundances. Notice endmembers. The matrix A ∈ ℜ m× + n

that all elements of b, A, and x are constrained to be positive and

∑x

i

= 1 . Typically m>n so we are dealing with an

i =1

Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XI, edited by Sylvia S. Shen, Paul E. Lewis, Proceedings of SPIE Vol. 5806 (SPIE, Bellingham, WA, 2005) · 0277-786X/05/$15 · doi: 10.1117/12.605670

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

711

overconstrained linear system of equations. This is called in the literature [1] the linear mixing model (LMM). There are other potential models that are discussed in [1] and references therein.

3

ABUNDACE ESTIMATION

3.1

Literature Review

For most applications, the measurement noise w in (1) is assumed to be independent identically distributed white Gaussian noise. That is w~N(0,σ2I) 2 where I is the m×m identity matrix and σ is the noise variance. The maximum likelihood estimate of x based on b is then given by

xˆ = arg min Ax − b x

2

(2)

2 n

subject to x ≥ 0 and ∑ xi = 1 i =1

where xˆ is the estimate of x, and || ||2 is the Euclidean norm. We will refer to (2) as the fully constrained unmixing problem. We could have arrived at LS minimization by plain curve fitting with no statistical considerations. This deterministic formulation allows us to look at the problem of estimating x as a problem of minimizing a measure of the difference between the measured pixel and the prediction of the pixel value based on the selected endmembers. In the literature, solutions to the unmixing problems are presented that partially enforce the abundance constraints. If no constraints are enforced in (2), the unconstrained least squares problem results for which direct solution is available see [2]. (3) xˆ = ( A T A ) −1 A T b ULS

Its main disadvantage in this approach in not having the constraints negatives values for the abundance could be obtained, having no physical meaning. In addition, the resulting estimate may not satisfy the sum to one restriction. We refer to this method as the Unconstrained Least Square method (ULS). A direct solution is also available for the sum-to-one constrained problem (see also [2]). Given by (2):

xˆ STO = xˆ ULS + (A T A )−1 λ

(4)

1 - 1x where λ = T T ULS−1 1 (A A) 1

[

]T

ˆ ULS is the solution of the unconstraint problem (2) and 1 = 1 1 ... 1 is a vector of ones of dimension n. where, x This approach still very simple to implement but the resulting abundance estimates may be negatives resulting in a solution with no physical meaning. For reference purpose, we refer to this method as the Sum to One Least Square (STOLS). In the non-negative least squares (NNLS) problem, only the positive constraint is enforced. For the NNLS problem only iterative methods can be used. The most used algorithm for the NNLS is an active set strategy of Lawson and Hanson described also in [2]. Multiplicative algorithms for the solution of the positively constrained abundance estimation problem are presented in [3]. To solve the fully constrained problem, [4] proposes to compute the abundance estimate as a Least Square with a penalty approach. We refer to this algorithm as Chang Fully Constraints (CFC) and refer to [4] for more details.

712

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

3.2 Abundance Estimation as a Least Distance Minimization Problem Our approach solves the fully constrained problem (2). To derive the proposed algorithm, we will transform (2) to a least distance minimization problem and use a Theorem presented in [2] to solve it using duality theory. The first transformation is based on the QR decomposition of the endmember matrix A given by

R  A = Q  0 where Q is an orthonormal matrix and R is an nxn upper triangular matrix. Define

 c1  QT b =   c 2  where c1 is an nx1 vector. Using this transformation, (2) can be transformed to [2]

xˆ = arg min Rx − c1 x

2

(5)

2 n

subject to x ≥ 0 and ∑ xi ≤ 1 i =1

We relaxed the sum to one to sum less than or equal to one since it allows the inclusion of a dark endmember easily. It is easily to shown, see [5], that (2) can be transformed to an inequality constrained similar to (5). Define (6) z = Rx − c 1

Substituting (6) into (5), we get

zˆ = arg min z z

2

(7)

2

subject to Gz ≥ g where

I 0  I  G =  T  R −1 , g =   −  T  R −1c1 1  1  1  This is a least distance (LD) problem [2]. The solution of the least distance problem using a dual method is based on the following equivalence between a LD and a particular NNLS problem derived in [2]. Theorem: Consider the LD problem (7). Let u be the solution of the NNLS problem

min Eu − f u

2

(8)

2

subject to u ≥ 0 where

G T   0 E =  T , f =   1 g  Let the residual corresponding to the solution be

r = Eu − f

(9)

and put σ=||r||2. If σ=0, then the constraints g≤Gz are inconsistent and (7) has no solution. If σ≠0, then the vector defined by

zˆ = [zˆ 1 L zˆ n ], zˆ j = − r j / rn+1 ,

j = 1,..., n

is the unique solution to (7). Proof: See Lawson and Hanson [2].

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

713

Once z is obtained, the abundance estimates are obtained from (6) as follows

xˆ = R −1 (zˆ + c1 ) = xˆ ULS + R −1zˆ

(8)

To solve (9), we can use any of the algorithms for NNLS described in [2] and [3]. The following pseudo-code describes our algorithm: Non Negative Sum Less or Equal to One (NNSLO) Algorithm:

R  0

 c1  . c 2 

1.

Compute the QR decomposition of A = Q   , compute c = Q b = 

2.

ˆ ULS = R c1 . Compute x

T

-1

T

[

T  T T 3. Let G = [I 1] , f = [0 1] , g = − x ˆ ULS −1  , E = G T  ˆ ULS 1 x 4. Set r = Eu − f , where u is the solution of the NNLS problem (6) for E and f. 5. Let zi = − ri / rn +1 ,1 ≤ i ≤ n T

6.

4

T

gT

]

T

.

−1

ˆ = xˆ ULS + R zˆ . The solution is: x

EXPERIMENTAL RESULTS

To compare the abundance estimates using the algorithms previously described, we use a Hyperion Image taken over Enrique Reef in La Parguera, southwest Puerto Rico. The Enrique Reef data used for this experiment consists of first 92 bands from 435-890nm with a spatial resolution of 30 meters. The data was pre processed with Fast Line-of-sight Atmospheric Analysis of Spectral Hypercubes (FLAASH) to remove atmospheric effects in the data. Figure1 shows an aerial photo of La Parguera and highlights the location of the reef.

Figure 1 Aerial Photographic of La Parguera, Lajas, Puerto Rico

Figure 2 shows a color composite image using the Enrique Reef HSI data and Figure 3 shows an image taken with the IKONOS sensor at 1 meter resolution. Clearly, we can see the degradation caused by the low spatial resolution of the Hyperion sensor (Figure 2) when compared with the IKONOS image (Figure 3). In addition, Figure 4 shows the spectral response of the selected endmembers: thalassia (sea grass), reef flat, rhizophora (mangrove), sand lagoon and sea water. The Pixel Purity Index (PPI) method of ENVI was used to extract the endmembers from the Enrique Reef image.

714

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

Figure 2 Enrique Reef Image from Hyperion Sensor

Sand

Sea Water

Sea Grass

Reef Crest Mangrove

Figure 3 Enrique Reef Image taken by IKONOS Sensor (1-meter resolution)

Figure 4 Enrique Reef Endmember

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

715

4.1 Abundance Estimation The algorithms ULS, STOLS, NNLS, CFC and NNSLO were implemented in ENVI and tested with the Enrique Reef data. Figure 5 shows a classification map of La Parguera area done by the Biogeography Program of the National Oceanic and Atmospheric Administration (NOAA) National Centers for Coastal Ocean Science, and is part of the Center for Coastal Monitoring and Assessment (CCMA) [7]. The method for producing the classifications maps is using aerial photographs to determine what kind of habitats, objects or structures are in the images, we refer to [7] for more information about the classification process.

La Parguera

Figure 5 La Parguera Classification Map

Figure 6 shows the abundance estimation results obtained for sea grass with ULS, STOLS, NNLS and NNSLO. The region estimated for sea grass is similar to that visually identified in the IKONOS image. We can seenegative values were estimates obtained by the ULS and STOLS algorithms. Abundance estimation results for sand with the ULS, STOLS, NNLS and NNSLO are shown in Figure 7. We also obtain negative abundances estimates with ULS and STOLS. Also, we can observe that the results obtained with NNSLO are closer to the IKONOS image. The results of the abundance estimates for sea water, mangrove and reef head endmembers were similar behavior with all algorithms. The NNSLO estimates are closer to the IKONOS image.

716

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

5

5

0

0

5

5

0

0

5

5

5

10

15

20

25

(a)

30

35

40

45

5

10

15

20

(c)

25

(b)

30

35

40

45

40

45

(d)

>100%

100-81% 40-21%

80-61% 20-0%

60-41% negative

Figure 6 Sea Grass Abundance Estimates with (a) ULS (b) STOLS (c) NNLS (d) NNSLO 5

0

5

0

5

5

10

15

20

25

(a)

30

35

40

5

45

10

(c)

15

20

25

(b)

30

35

(d)

>100%

100-81% 40-21%

80-61% 20-0%

60-41% negative

Figure 7 Sand Abundance Estimates with (a) ULS (b) STOLS (c) NNLS (d) NNSLO

Figure 8 shows a comparison of the abundance estimates of CFC and NNSLO for the different endmembers. We can observe that both algorithms satisfy the positivity and the sum to one constraint. The estimates of sea grass, reef crest, mangrove and sea water results for both algorithms are very similar.

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

717

CFC

NNSLO

Sand

Sea Grass

Reef Crest

Mangrove

Sea Water

100-81%

80-61%

60-41%

40-21%

20-0%

Figure 8 Comparisons in Abundances Estimates between CFC and NNSLO

The sand abundance estimates obtained with the NNSLO are closer to the IKONOS than CFC. All comparisons are visual.

718

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

5

CONCLUSIONS

This paper presented an algorithm to solve the fully constrained abundance estimation problem. The algorithm is based on transforming the fully constrained abundance estimation problem to an equivalent least distance problem and using a Theorem from Lawson and Hanson that solves this problem using dual theory as a nonnegative least squares problem. We also present result of comparison of different algorithms: ULS, STOLS, NNLS and CFC in unmixing a Hyperion image of the La Parguera, Puerto Rico and compared them with results obtained with NNSLO. As the results shows, the NNSLO algorithm abundance estimates that gives the best visual agreement with a high spatial resolution IKONOS image. ACKNOWLEDGEMENTS This work was sponsored primarily by the NASA University Research Centers Program under grant NCC5-518. The research performed was partially funded by the Center for Subsurface Sensing and Imaging Systems at the University of Puerto Rico-Mayaguez sponsored by the Engineering Research Centers Program of the US National Science Foundation under grant EEC-9986821. In addition, w like to thank to Francisco Rivera-Maldonado and Dr. Raul Torres for providing the IKONOS image with the segmentation. REFERENCES [1] [2] [3]

[4] [5] [6] [7]

N. Keshava and J.F. Mustard, “Spectral Unmixing.” In IEEE Signal Processing Magazine, pp. 44-57, January 2002. C.L. Lawson and R.J. Hanson, Solving Least Squares Problems, Prentice-Hall, 1974. M. Velez-Reyes, S. Rosario, A. Puetz, R. B. Lockwood, “Iterative algorithms for unmixing of hyperspectral imagery.” In Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery IX, Proceedings of SPIE Vol. 5093, April 2003. Chang C. I.; Ren H.; et. al. ‘‘Subpixel Target Size Estimation for Remotely Sensed Imagery". In Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery IX, Proceedings of SPIE Vol. 5093, April 2003. Rosario, S. Iterative Algorithms for Abundance Estimation on Unmixing of Hyperspectral Imagery. Master of Science Thesis, University of Puerto Rico at Mayagüez, May 2004. M Vélez-Reyes and S. Rosario, “Solving Adundance Estimation in Hyperspectral Unmixing as a Least Distance Problem” In Proceedings of International Geoscience and Remote Sensing Symposium, September 2004. Kendall, M.S., M.E. Monaco, K.R. Buja, J.D. Christensen, C.R. Kruer, and M. Finkbeiner, R.A. Warner. 2001. “Methods Used to Map the Benthic Habitats of Puerto Rico and the U.S. Virgin Islands”, http://biogeo.nos.noaa.gov/projects/mapping/caribbean/startup.htm

Proc. of SPIE Vol. 5806

Downloaded From: http://proceedings.spiedigitallibrary.org/ on 01/03/2016 Terms of Use: http://spiedigitallibrary.org/ss/TermsOfUse.aspx

719