Blind Source Separation Applied to Spectral Unmixing - Muscle

0 downloads 0 Views 591KB Size Report
Blind Source Separation Applied to Spectral. Unmixing: ... we analyze the performance of some of these measures on an experimental basis. It is shown that one ...
Blind Source Separation Applied to Spectral Unmixing: Comparing Different Measures of Nongaussianity Cesar F. Caiafa1 , Emanuele Salerno2 , and Araceli N. Proto1,3 1 Laboratorio de Sistemas Complejos. Facultad de Ingenier´ıa - UBA, Av. Paseo Col´ on 850. 4to. piso, Ala sur (1063) Capital Federal, Argentina 2 Istituto di Scienza e Tecnologie dell’Informazione - CNR Via Moruzzi, 1, I-56124 Pisa, Italy 3 Comisi´ on de Investigaciones Cient´ıficas de la Prov. de Buenos Aires, Av. Paseo Col´ on 850. 4to. piso, Ala sur (1063) Capital Federal, Argentina [email protected],[email protected],[email protected]

Abstract. We report some of our results of a particular blind source separation technique applied to spectral unmixing of remote-sensed hyperspectral images. Different nongaussianity measures are introduced in the learning procedure, and the results are compared to assess their relative efficiencies, with respect to both the output signal-to-interference ratio and the overall computational complexity. This study has been conducted on both simulated and real data sets, and the first results show that skewness is a powerful and unexpensive tool to extract the typical sources that characterize remote-sensed images. Keywords: Blind Spectral Unmixing, Dependent Component Analysis, Hyperspectral Images, Maximum non-Gaussianity, Unsupervised Classification.

1

Introduction

Airborne and satellite-borne hyperspectral sensors have given new possibilities for remote sensing applications. For example, the availability of accurate and highly resolved spectral information has enabled remote geological surveys of large areas. Other important possibilities have been provided by the very high spatial resolution often available. Unfortunately, spatial and spectral resolutions are intrinsically conflicting requirements, since the total energy that a sensor can collect from a single spatial resolution cell is directly related to the cell area, the sensor bandwidth, and the integration time. Then, it is apparent that, for a fixed integration time and a fixed (narrow) band, the total energy coming from a very small cell can be insufficient to ensure the minimum required signal-to-noise ratio. In the case of moving sensors, the integration time must satisfy precise requirements, and cannot be increased at our will; thus, the spatial resolution can be much poorer than optically attainable, and even spectrally distinct materials B. Apolloni et al. (Eds.): KES 2007/ WIRN 2007, Part III, LNAI 4694, pp. 1–8, 2007. c Springer-Verlag Berlin Heidelberg 2007 

2

C.F. Caiafa, E. Salerno, and A.N. Proto

become undistinguishable when contained in a single resolution cell. Spectral unmixing consists in evaluating the percent occupations (or relative abundances) of the individual materials (or endmembers) within a resolution cell on the basis of their reflectance spectra. This problem would not be so complicated if all the endmember spectra were known but, unfortunately, this is a very unlikely situation. Extracting endmember relative abundances with no knowledge about their spectra is called blind spectral unmixing, and has been addressed recently by a number of authors [1,2,3]. One of the most popular attempted approaches is based on independent component analysis and assumes that the endmember abundances can be modeled as independent random variables. This assumption is strong, and not realistic for many reasons [1,3]. This is why we proposed a new dependent component analysis strategy for blind spectral unmixing (MaxNG), which is based on the local maximization of nongaussianity without any independence constraint [4]. In the particular appplication of spectral unmixing, furthermore, blind separation does not suffer the scale ambiguity that is common to all the results of blind processing in the absence of additional constraints [5,6]. Our first experimental results [4,6] have shown that, applied to spectral unmixing, our approach is more effective than some classical independent component analysis techniques. There is a computational complexity issue that needs to be addressed, for which we propose a solution in [6]. Another aspect to be assessed is the efficiency of the nongaussianity measure we adopt in performing our search of the optimal endmember separation. MaxNG is based on the Euclidean distance in the L2 space. Other measures could be chosen, such as entropy, negentropy, skewness, kurtosis, higher order moments, etc. In this paper we analyze the performance of some of these measures on an experimental basis. It is shown that one of the simplest measures, skewness, gives excellent results with a very low computational cost. In Section 2, we briefly describe MaxNG as applied to blind spectral unmixing. In Section 3, we summarize the different nongaussianity measures we are taking into account and, in Section 4, we show our experimental results. A brief discussion closes the paper.

2

Blind Spectral Unmixing and MaxNG

Let us assume we have a hyperspectral reflectance data stack x(i) that, for each pixel i, is constituted by a real M -vector x, ideally given by the following linear mixture model: x(i) = As(i) + n(i),

i = 0, . . . N − 1

(1)

where s(i) is a real P -vector containing the fractional endmember abundances for pixel i, A is a pixel-invariant M × P matrix whose (k, l)-th entry represents the value of the l-th endmember spectrum at the wavelength of the k-th channel, and n(i) is a real M -vector of additive noise (independent of s(i)) at pixel i. Matrix A

Blind Source Separation Applied to Spectral Unmixing

3

is assumed unknown, and the elements of s are considered as dependent random variables subject to the condition P −1 

i = 0, . . . N − 1

sl (i) = 1,

(2)

l=0

This is because, in this model, the spectrum of the light scattered by each pixel is obtained by the sum of the individual endmember spectra, each weighted by the percent occupation of the related endmember. Spectral unmixing consists in estimating vector s at each pixel from the available data stack x. In the present case, we do not assume any knowledge of endmember reflectance spectra, that is, of matrix A, and we talk of blind spectral unmixing. The idea under MaxNG [4] is that a linear blind spectral unmixing can be achieved by combining the data through an M -vector d and finding the local maxima of a measure of nongaussianity for the transformed variable (i), where x  is a whitened version of the data stack (see also [7]). z(d; i) = dT x In hyperspectral imaging, the number of available channels M is usually much larger than the number of endmembers. For this reason, the size of the whitened  does not need to be the same as the one of vector x; for example, vector x when P endmembers are present in the scene, Eq. (2) implies that the effective data dimensionality is at most Q = P − 1. If we find a vector d∗ such that the nongaussianity of variable z(d∗ ; i) is a local maximum, then z(d∗ , i) in all the N pixels is a scaled version of one of the endmember abundance maps. The performance of MaxNG will depend on the particular nongaussinity measure adopted. An overview of possible options is reported in the next section. Under constraint (2), we have shown that it is also possible to estimate the scaling parameters, thus accomplishing unambiguous spectral unmixing [6].

3

Non-gaussianity Measures

Given a zero-mean and unit-variance random variable z, we can devise many ways to measure its nongaussianity. We report here the measures we adopted to test the performance of MaxNG. – Euclidean distance in L2 [4]:  ΓN G =

+∞ −∞

[fz (z) − Φ(z)]2 dz

(3)

where fz (z) is the pdf of z and Φ(z) is the standard Gaussian. This distance has already been used for blind separation of dependent sources [4,6]. – Negentropy [7]:  ΓN egentropy = HGauss −

+∞

−∞

fz (z) log(fz (z))dz

(4)

4

C.F. Caiafa, E. Salerno, and A.N. Proto

√  where HGauss = log 2πe is the Shannon differential entropy associated to the standard Gaussian. This measure, and some approximations thereof, have been exploited for independent component analysis. – Central moments of order P > 2 [7,8]. In particular, for P = 3, we have the skewness, whose magnitude is a measure of asymmetry for a pdf. A measure of nongaussianity could be the squared skewness:  +∞ 2 2 3 z fz (z)dz (5) Γm3 = m3 = −∞

Among the other higher-order moments and cumulants, kurtosis has already been used successfully for independent component analysis [7]. We also tried with the sixth-order moment:  +∞ Γm6 = m6 = z 6 fz (z)dz (6) −∞

Using the Parzen windows nonparametric pdf estimation technique [9], measures (3) and (4) can be estimated from a set of N samples [11,4]. The drawback is  that the related computational complexity is O N 2 . FFT-based algorithms [10] can reduce the complexity to O (N + NΔ log NΔ ), with NΔ < N . Conversely, higher order moments show the great advantage of being easily computable by averaging, with complexity O (N ): N −1 1  P z (i) mP ≈ N i=0

(7)

To maximize the distance functions, it is necessary to differentiate them with respect to the search parameter d. The moment-based measures can easily be differentiated with complexity O (N ): ∇d mP ≈

N −1 P  P −1 z (i) x(i) N i=0

(8)

On the other hand, the L2 and the negentropy distances are more complicated to differentiate. Again, the related complexity can be reduced by FFT-based techniques [6].

4

Experimental Results

Experiments with Synthetic Data: To synthesize P source variables with the features summarized in Section 2, we first sampled P independent random P −1 variables w0 , w1 ,.., wP −1 from a common pdf, and then let sk = wk / i=0 wi , 1 with k = 0, . . . P − 1. 1

Another possibility for generating sources is to use the Dirichlet distribution (Multivariate Beta) as done in [3].

Blind Source Separation Applied to Spectral Unmixing

5

We tested the different nongaussianity measures for several SNRs2 and numbers of samples N for a set of P = 3 synthetically generated sources. A total number of channels M = 102 and randomly generated mixing matrices A were used. Mixtures xk were generated according to the linear model (1), and uncorrelated Gaussian noise with channel-independent SNR was introduced.  For the different cases, the nongaussianity measures of variable z(d) = dT x  being the 2D whitened mixture vector3 . Therefore, d were calculated, with x can be written in terms of a single parameter θ (see [4]): dT = [cos θ, sin θ]. To illustrate our results, in Fig. 1, we report the nongausianity measures as functions of θ. Each plot is completed with vertical lines placed at the values of θ for which the normalized sources are recovered from the 2D whitened data . The figure shows four cases with different values of SNR and N , and vector x different nongaussianity measures. It is worth mentioning that the measure based on kurtosis, although it proved to be very efficient in independent component analysis, did not give good results in this application and was not reported in our plots. Possible theoretical reasons for this result will be studied in the future. We observe that, for high SNRs, all the local maxima are very close to their expected positions. Additionally, all the considered measures showed some robustness to noise. As could be expected, the most sensitive to noise seems to be Γm6 . As an index of the quality of separation for each estimated source sˆl , we assumed the usual signal-to-interference ratio (SIR), defined as SIRl = 10 log10 (var(sl )/var(ˆ sl − sl ))

(9)

In Fig. 2, a comparison among the results of applying different measures is shown. For each measure, we averaged the results over a large set of simulations: the Mean SIR values (Fig. 2.a) were calculated for an SNR range −20dB to 60dB over a set of 60 source estimations for each SNR value. As a reference, we compare these results with the case of using the optimal linear estimator D (DA = I). In this case, the estimation error is only caused by the additive Gaussian noise. To analyze how noise affects the local maximization, we defined the efficiency of separation as the ratio of the number of sources successfully detected over the total number of cases. We assume that a source has been recovered correctly if a minimum SIR of 8 dB is obtained when the noiseless vector As(i) is pre-multiplied by the obtained separating matrix D, i.e., ˆs(i) = DAs(i). The resulting efficiency values are shown in Fig. 2.b. Note that the negentropy and the sixth-order moment measures are more sensitive to noise than the L2 distance and the skewness measures. Experiments with Real Hyperspectral Images: We applied MaxNG with the L2 distance and skewness measures to analyze a piece of 75 × 75 pixels (N = 5625) taken from a real, radiometrically corrected, hyperspectral image 2

3



SN R(dB) = 20 log σσnx , where σx and σn are the standard deviations of the useful signal and the added noise, respectively. Note that the constraint (2) reduces the dimension of the source space to 2D.

6

C.F. Caiafa, E. Salerno, and A.N. Proto

b) N=512 and SNR=0dB

a) N=512 and SNR=50dB GNG GNegentropy

0.8

Gm3 Gm6

0.6

GNG

1.0

Normalized Index

Normalized Index

1.0

Detected Max

0.4 0.2

GNegentropy

0.8

Gm3 Gm6

0.6

Detected Max

0.4 0.2

0.0

0.0 0

20

40

60

80

100

120

140

160

180

0

20

40

60

80

q

c) N=4096 and SNR=50dB GNegentropy

0.8

Gm3 Gm6

0.6

120

140

160

180

d) N=4096 and SNR=0dB 1.0

GNG

Normalized Index

Normalized Index

1.0

100

q

Detected Max

0.4 0.2

GNG GNegentropy

0.8

Gm3 Gm6

0.6

Detected Max

0.4 0.2

0.0

0.0 0

20

40

60

80

100

q

120

140

160

180

0

20

40

60

80

100

120

140

160

180

q

Fig. 1. Nongaussianity measures based on L2 distance, negentropy, m3 , and m6 versus angle θ. All values were normalized in the range [0,1]. Vertical lines denote the expected positions of the local maxima.

Fig. 2. MeanSIR (a) and Efficiency (b) for NG measures based on L2 distance (NG), Negentropy, m3 and m6

from an urban area in Rome, Italy. It contains 102 data channels with center wavelengths ranging from 0.43μm to 12.70μm. These data were also accompanied by a pixel classification. Four supplied classes (grit, vegetation, bricks, roads, from left to right) are shown in the top row of Fig. 3. This classification has been

Blind Source Separation Applied to Spectral Unmixing

7

Fig. 3. Top, left to right: grit, vegetation, brick, road classes associated to our test 102band hyperspectral image. Middle: four fractional endmember abundances extracted by MaxNG with L2 -distance nongaussianity measure. Bottom: four fractional endmember abundances extracted by MaxNG with squared-skewness nongaussianity measure.

obtained by standard, nonblind, techniques, and assigns a class to each pixel. As explained above, spectral unmixing assigns each pixel with a set of fractional abundances. Thus, the binary classes supplied with the image data cannot be used as the proper ground truth. We only used them to visually appreciate the quality of our estimates. The center and bottom rows in Fig. 3 show our results obtained by Euclidean distance in L2 and skewness, respectively. It is apparent that the skewness-based results are very close to the ones obtained by the Euclidean distance, with the advantage of a shorter running time.

5

Discussion

This paper reports our first results in comparing different measures of nongaussianity as applied to our previously proposed MaxNG blind spectral unmixing technique for hyperspectral images. The Euclidean distance in L2 assumed as nongaussianity measure in our original formulation proved to be very efficient and insensitive to noise, but entails a considerable computational cost, especially when the data sets are very rich. Conversely, the measures based on higher-order moments or cumulants are characterized by a low computational cost (O (N )). This was our motivation for investigating their performance when used with MaxNG in the particular task of spectral unixing. The first unexpected result was that the kurtosis-based measure is uneffective in this application. A theoretical insight will be needed to explain this result. The skewness-based measure performed very well, with the mentioned computational advantage over the Euclidean distance. This result can be perhaps explained by observing that the endmember distributions

8

C.F. Caiafa, E. Salerno, and A.N. Proto

are normally very skewed, and this feature is better captured, of course, by a skewness-based measure. Our future research on MaxNG will include a more extensive experimentation of different nongaussianity measures, as well as a further investigation of the theoretical issues raised in [6]. In particular, we will try to relax some of the simplifying assumptions underlying the present data model. Acknowledgements. This work has been partially supported by the EU Network of Excellence MUSCLE (FP6-507752). C. Caiafa acknowledges financial support from Facultad de Ingenieria, Universidad de Buenos Aires, Argentina (Beca Peruilh). The 102-band hyperspectral image is courtesy of the Airborne Laboratory for Environmental Research at IIA-CNR in Rome, Italy. The authors are indebted to Lorenza Fiumi for helpful discussions.

References 1. Keshava, N., Mustard, J.: Spectral unmixing, IEEE Signal Process. IEEE Signal Process. Mag. 19(1), 44–57 (2002) 2. Berman, M., Kiiveri, H., Lagerstrom, R., Ernst, A., Dunne, R., Huntington, J.: ICE: A statistical approach to identifying endmembers. IEEE Transactions on Geoscience and Remote Sensing 42, 2085–2095 (2004) 3. Nascimento, J.M.P., Bioucas Dias, J.M.: Does Independent Component Analysis Play a Role in Unmixing Hyperspectral Data? IEEE Transactions on Geoscience and Remote Sensing 43, 175–187 (2005) 4. Caiafa, C.F., Proto, A.N.: Separation of statistically dependent sources using an L2 -distance non-Gaussianity measure. Signal Processing 86, 3404–3420 (2006) 5. Caiafa, C.F., Salerno, E., Proto, A.N., Fiumi, L.: Dependent Component Analysis as a Tool for Blind Spectral Unmixing of Remote Sensed Images. In: Proc. EUSIPCO 2006, Florence, Italy, September 4-8, 2006, pp. 4–8 (2006) 6. Caiafa, C.F., Salerno, E., Proto, A.N., Fiumi, L.: Blind Spectral Unmixing by Local Maximization of Non-Gaussianity. Signal Processing (submitted) (2006) 7. Hyv¨ arinen, A., Karhunen, J., Oja, E.: Independent Component Analysis. J. Wiley & Sons, New York (2001) 8. Cruces, S., Cichocki, A., Amari, S.-I.: The Minimum Entropy and Cumulants Based Contrast Functions for Blind Source Extraction. In: Mira, J.M., Prieto, A.G. (eds.) IWANN 2001. LNCS, vol. 2085, pp. 13–15. Springer, Heidelberg (2001) 9. Parzen, E.: On estimation of a probability density function and mode. Annals of Mathematical Statistics 33, 1065–1076 (1962) 10. Silverman, B.W.: Density Estimation for Statistics and Data Analysis. Chapman and Hall, New York (1985) 11. Erdogmus, D., Hild, K.E., Principe, J.C., Lazaro, M., Santamaria, I.: Adaptive blind deconvolution of linear channels using Renyi’s entropy with Parzen window estimation. IEEE Transactions on Signal Processing 52, 1489–1498 (2004)