Unsupervised linear spectral mixture analysis with AVIRIS data

Journal of Harbin Instiule of Technology (New Series), Vol. 12, No. 5, 2005Unsupervised linear spectral mixture analysis with A VIRIS dataGU Yan-feng' , YANG Dong-yun2 , ZHANG Ye'谷延锋,杨冬云,张晔(1. Dept. of Informnation Engineering, Harbin Institute of Technology, Harbin 150001，China;2. Dept. of Electronic Engineering, Heiongjiang Institute of Technology，Harbin 150050, China)Abstract: A new algorithm for unsupervised hyperspectral data unmixing is investigated ，which includes a mod-ifed minimum noise fraction ( MNF) tranformation and independent component analysis ( ICA). The modifiedMNF transformation is used to reduce noise and remove correlation between neighboring bands. Then the ICA isapplied to unmix hyperspectral images , and independent endmembers are obtained from unmixed images by using post-processing which includes image segmentation based on statistical histograms and morphological opera-tions. The experimental results demonstrate that this algorithm can identify endmembers resident in mixed pix-els. Meanwhile, the results show the high computational fficiency of the modified MNF transformation. Thtime consumned by the modified method is almost one fith of the traditional MNF transformation.Key words: spectral mixture analysis; minimum noise fraction; independent component analysis; linear mix-ture model; adaptive subspace decompositionCLC number: TN911. 73Document code: AArticle ID: 1005-91 13(2005 )>05 0471-06The recent development of hyperspectral sensorsfractions in a mixed pixel are called fractional abun-has shown potential for monitoring the earth’s surface,dance 3。Up to now, two models have been proposedconducting resource surveys and target detection. Hy-to describe spectral mixture. One is the macro-spectralperspectral sensors can image an area in hundreds ofmixture that models a mixed pixel as a linear combina-different continuous wavelength ranges. An example oftion of signatures in the pixel with relative concentra-such a hyperspectral sensor is the airborne visible/in-tions. The other model is called a microscopic or inti-frared imaging spectrometerAVIRIS ) with 224mate mixture and is a nonlinear mixing model of signa-bands' . Thus hyperspectral data provide a dense sam-tures within the pixel. Nevertheless, the second modelpling of the spectral signatures of land covers, allowingcan be linearized by a semi-emprical method 5。Thea better discrimination among similar classes than tradi-Linear Mixture Model ( LMM) is built on such an as-tional multispectral scanners2 .sumption that mixed pixels are linear combinations ofIn hyperspectral images, mixed pixels are a mix-several distinct materials. Although it cannot describeture of more than one distinct material, and they existall practical cases of the spectral mixture, LMM hasfor the following reasons 3. Firsl, the spatial resolu-shown better performance in many applications. How-tion of the sensor limits the physical size of the imageever, a main disadvantage in most known linear unmix-pixel and makes a single pixel tend to include differenting algorithms is the need for a priori knowledge of thematerials so that the resulting spectral measurement willmaterial signatures present. This is very difficult to ac-be composites of individual spectra. Second, mixedquire in practice.pixels can occur when distinct materials are combinedTraditional MNF transformation is a linear transfor-into a homogeneous mixture which is independent of themation that maximizes the signal-to-noise ratio (SNR)，spatial resolution of the sensor. Third, energy transfer thus arranging images in terms of decreasing imageinto and away from neighboring instantaneous field ofquality. Like principal components analysis ( PCA)，views ( IFOVs) leads to further pixel ambiguity [4)the MNF transformation can be used to reduce noiseTo solve the problem of mixed pixels, it is neces-and remove correlation between neighboring bands.sary to investigate spectral mixture analysis ( SMA).When the MNF transformation is carried out in theThe basic premise of modeling mixture is that within awhole space of hyperspectral images, the MNF transfor-given scene, the surface is dominated by a small num-mation. requires a huge amount of computation. Mean-ber of distinct materials, which have relatively constantwhile ,中国煤化工n the neighboringspectral properties and are called endmembers. Thebands_Inds far from eachFYHCNMHGReceived 2003 -06 -26.Sponsored by the National Natural Science Foundation of China( Grant No.60272073)..471●Journal of Harbin Instiule of Technology (New Series), Vol. 12, No. 5, 2005other, the useful transformed features provided by thetraditional MNF，which performs in the entire data1 Modified MNF Transformationspace and ignores local correlation, are insufficient forthe following processing. Therefore, it is very necessa-Traditional MNF transformations are realized inry to modify the traditional MNF transformation in orderthe whole original data space, and often require hugeto reduce the amount of data to be processed and to im-amounts of computation.In addition, the direct pro-prove computational efficiency. .cessing in the original data space also means that theICA is a technique that comes from blind sourceuseful transformed features are insuficient, due to aseparation in the signal processing community, and hasweaker correlation of bands in terms of the whole datarecently gained attention in speech enhancement, tele-space rather than that of local neighboring bands. Incommunications，medical signal processing, etc. Withtermns of these factors, a modified method is provided torespect to image processing, Bayliss etal (81 and Tu etimprove the quality of the MNF transformation and itsal ' ”applied ICA to mixture analysis in hyperspectralcomputational efficiency by introducing the spectralimages. Furthermore, Lee et al 240 investigated unsu-subspace-processing scheme. The following describespervised image classification， segmentation and en-the subspace-processing scheme and the traditionalhancement using ICA mixture Models. By assumingMNF transformation and introduces the modified meth-statistical independence between different sources， 8linear non-orthogonal coordinate system is found by1.1 Spectral Subspace-processing SchemeICA in multivariate data determined by second-andFirstly，the spectral-dimensionality feature vectorhigher-order statistics. The goal of ICA is to linearlyis built by using a band of hyperspectral images as atransform the data so that the transformed variables arerow vector. Then the correlation matrix of the featureas statistically independent from each other as possi-vector is calculated, and the visualization of the corre-ble. Therefore, when it is used to SMA, ICA can real-lation matrix is shown in Fig. 1. From this the blockingize unsupervised linear spectral unmixing in blind cir-feature can be found, and it clearly demonstrates thecumstances without priori knowledge. Thus it can over-strong correlation between neighboring bands. There-come the drawback of LMM.fore, it is more suitable to process the hyperspectralIn this paper, a new unsupervised linear spectraldata in subspace with stronger correlation. An adaptivemixture analysis algorithm is proposed, which can dis-subspace decomposition ( ASD) scheme is applied tocriminate，identify and quantify independent endmem-reduce the dimensionality of hyperspectral datalo. ' andbers in a mixed pixel without any priori knowledge.is an efficient approach to reduce the huge arnount ofThe algorithm consists of lwo main procedures: a modi-data for processing by the MNF transformation. Thefied minimum noise fraction ( MNF ) transformation andASD method is adapted in the modified MNF transfor-independent component analysis ( ICA). By introdu-mation.cing the subspace-processing scheme， the modifiedMNF transformation is carried out in the subspace andused to reduce noise and decorrelate data. Then ICA isused to model the mixture in hyperspectral images ，separate individual materials in mixed pixels in an un-supervised manner and determine their mixing propor-tion. After the ICA mixture model is resolvedby usingthe infomax algorithm" ， threshold segmentationbased on the statistical histogram and morphologicaloperations are used to find endmembers.The paper is organized as follows. In Section 1,the improvement method for MNF transformation is pro-Fig.13D Correlation map with blocking feature betweenposed by using spectral a subspace-processing scheme ,neighboring bands of hyperspectral imagesafter introducing the traditional MNF transformation. InSection 2, the ICA algorithm is introduced for SMA,According to correlation between neighboringand it is performed on the transformed data with the bandsa n gloha theshold narameter is delemined,modified MNF without priori knowledge. In Section 3,base中国煤化工data space P witha framework of byperspectral image unmixing by usingN-dCNMH Gto I subspacesP,(j =the modified MNF transform and ICA is provided. TheYHsionaliydescription of hyperspectral image data and numericalexperiments are given in Section 4. Finally, we con-P=∪P,clude this paper in Section 5 with a short summary.Notice that different subspace has different num-●472.Journal of Harbin Instiute of Technology (New Series), Vol. 12, No. 5, 2005bers of dimensionality and that similar spectral charac-pixels contain a similar signal but also random isolatedteristics exist in the bands of each subspace.noise. A simple and rapid method, called shift differ-1.2 Minimum Noise Fraction Transformationence, is applicable for estimating E。，which can beTraditional MNF transformation was first proposeddescribed asby Green et al (12]. Gordon [1] generalized the MNFN(x) = D(x) -D(x +δ),transformation by introducing powers of each band asS(x) = S(x +δ),new bands before performing the MNF transformation.where δ is an appropriately determined step length.To eficiently compute the MNF transformation matrix,Here, it is necessary to note that the best noise estima-Roger 14) investigated a matrix decomposition method,tion is obtained by using the shift difference statisticswhich simplifies the computation of the MNF transfor-from a local homogeneous area rather than from themation to some extent.Generally, the MNF transformation includes twowhole image.After obtaining the noise covariance 2 y，we ap-steps. The first step is to remove band-to-band correla-tions and rescale the noise data with unit variance.ply matrix eigenvalue decomposition with 2、to com-This step requires that the noise covariance matrix ofpute its ortho-normalized eigenvector matrix Hx, andthe data be known or estimated.Consider an N-bands hyperspectral data set withthe diagonal matrix of its eigenvalues, E. Thus, wegray levels D(x) = {D、(x),.,Dv(x)IT, where xcan obtain noise-whitening matrix WN, which can res-gives the coordinate of the sample. The MNF transfor-cale the noise data to have unit variance, namely,mation can be described by the following expressionw、E、W、=I,(1)Z(x) = A'D(x)，where WN = H、ET , and I is identity matrix.The second step of the MNF transformation is awhere Z(x) is the transformed data with dimensionalitystandard PCA transformation of the normalized data.Mi.e. z(x) = {Z;(x)，2(x)， .,.(x)1,Combining the first step with the second step, weand A is the transformation matrix with size M X N.can obtain the final expression of linear transformationTo obtain the transformation matrix A, a process-matrix A as the followinging step similar to the following is adapted. First, theA 2、2。= AA,noise covariance statistic is estimated from the hyper-spectral data. Then a whitening matrix is constructed towhere A is a diagonal matrix of the eigenvaluesλ,(i =whiten the data by using the noise statistic. Finally, by. 1,2,.,N)of 2。Note that the second step of thecombining the result from the first step of processingprocedure is performed on the noise-whitened data, inwith the PCA transformation, the final MNF transfor-which the noise covariance is the identity matrix.mation matrix A can be obtained. The detailed processTherefore, the SNR is decided by the signal covarianceis described as follows.二。By ordering theλ; so thatλ≥λ2≥.≥λv，Assume that the noise model is additive i. e.D(x) = S(x) +N(x) ,the MNF transformed data are arranged in bands of in-where S(x) and N(x) are the uncorrelated signal andcreasing noise fraction. Thus the transformed data arenoise components. Thus, if the mean of D(x) is as-arranged in terms of decreasing SNR.sumed to be zero, the covariance matrix 三。of D(x)2 ICA for Unmixing Hyperspectral Imagemay be expressed bycov|D(x)1 = E,= 2。+ E,The goal of ICA is to linearly transform the data sothat the transformed variables are as statistically inde-whereE, and E、are the covariance matrices ofpendent from each other as possible. Generally, ICAS(x) and N(x)， respectively.The noise fractionrequires two basic premises. One is that the source(NF) of the i-th band is defined ascomponents are statistically independent. The otherNF = var{N,(x)}/var{D,(x)}.premise is that, at most, the probability density func-The MNF transformation converts observation datation of one independent source component has Gaussian.from the original coordinate system into the new sys-distribution58.hose premises can be satisfied intem, in which the transformed data is the linear combi-termns. The Informax algo-nation of the original data.rithm中国煤化工ach to ICA181.In order to obtain the transformation matrixA, theYHC NM H G) ' are n stisticallynoise covariance matrix > has to be known or esti-independent endmembers, and x = (x,x2.,x)mated. Assume that each pixel of the hyperspectral im-are n observations and are the linear combination of in-age contains both signal and noise, and that adjacentdependent endmembers s. Notice that x are the trans-●473●Journal of Harbin Instiute of Technology (New Series), Vol. 12, No. 5, 2005formed data obtained by the modified MNFT, practical-tained by learning with a neural network. Thus ，inde-ly. Thus, the relation of x and s is expressed bypendent endmembers are obtained by using formulax =Ms(2) after the separating matrix is obtained.where M is the mixing matrix with full rank. In the ICAalgorithm, the endmembers s and the mixing matrix M3 Unsupervised Linear Spectral Unmixingare unknown. Under these conditions and the two as-sumptions mentioned above, the approach to resolveAccording to the discussion above, it is found thatthe ICA mixture model is to find a linear mapping Wthe modified MNF transformnation can be used to rewith n dimensionality so that the unmixed endmembersmove noise, rescale noise with unit size and determiney satisfy the following expressionthe inherent dimensionality of hyperspectral data. They=Wx=WMs2)ICA can be used to unmix hyperspectral images withwhere components ofy, i. e. the estimated independentlinear models in blind circumstances. Therefore, weendmembers，are statistically independent, and thcombine both the modified MNF transformation and thematrix W is a separating matrix.ICA to construet an unsupervised linear spectral mix-In the ICA algorithm, maximizing the nongaussi-ture analysis algorithm. The system diagram is shownanity of the data allows the independent endmembers toin Fig.2.be obtained. The central limit theorem is the basis formaximizing the nongaussianity of the data. Generally ,AVIRIS datathere are two methods used to measure the nongaussi-anity of the data: kurtosis and infornation entropy.The drawback of kurtosis is that it is sensitive to outli-ers in the data. The idea of the informax algorithm is toSubspace PartitionStep 1maximize mutual information between s and y, in whicha crucial process is to find the appropriate separatingSubspace 1Subspace Imatrix in an iterative manners..According to formula (2) , the probability density↓function of the observation data x is described byp(x) =| det(W)I p(y),(3)hit differenceShift diferencewherep(y) = Ip(y) is defined as an estimate ofData W hiteningData Whiteningindependent endmembers’pdf. Logarithmic operationis performed on both sides of formula (3) and the like-lihood function L(y, W) is obtainedL(y,W) = logp(x) = log | det(WV) 1+PC/PCA2 logp(y,).Using the maximum likelihood estimate ( MLE)StepsICA misture modelwill maximize formula (4) and obtain the separatingmatrix. The stochastic gradient ascent-learning rule forW can be obtained by the following expressionPost-processing△W =aL(y,W)=oW[(W")-' -o(y)x" ,ifdet(W) > 0Endmenmbers1- (W")- -φ(y)x ,ifdet(W) <0ap(y)/dyap(y)/dy_Fig.2 System diagram for unsupervised linear spectral un-where p(y) =-= [-p(y)p(y)，-mixingap(y2)/ayap(y.)/2yThus, the leamingIn the unsupervised linear spectral mixture analy-equation with iterative form is obtained by means of thep(x2)p(yn)sis algorithm, the whole hyperspectral data space isfin中国煤化工aces according to thelearning rulecorrands. Then the shiftWk1=W。+aOWdifMHCN M H Glate noise dala covari-where a is the learning rate used to control the learningance from a homogeneous region, and whiten the noiserule and the convergence speed. When numerical iter-that exists in hyperspectral images by using the noiseation converges, the final separating matrix W is ob-●474●.Journal of Harbin Institute of Technology (New Series), Vol. 12, No.5, 2005covariance. In the following, standard PCA transforma-homogenous region in the data subspace ，which corre-tion is performed on the whitened data.sponded to a soybeans-min area in the ground and in-After the modified MNF transformation, we obtaincluded 580 pixels. After obtaining the estimated noisethe decorrelated hyperspectral images with high SNR,covariance, we carried out the data-whitening processof which dimensionality is reduced. Then we can un-with it. Then, the PCA transformation was performedmix the obtained data using the ICA mixture model. Inon the whitened data in each of the subspaces.order to resolve the mixing matrix and to obtain unmix-In order to verify the computational efficiency ofing materials with the ICA mixture model, we adopt thethe modified MNF transformation， we compared theinfomax algorithm. In this way, the unmixing materialstime consumed by both the traditional and the modified:an be obtained by inverse mapping with the mixingAINF transformation. The time consumed by the tradi-matrix. In order to reduce the disadvantageous effect oftional MNF transformation was more than 102 seconds.the residual , we segment the obtained unmixing materi-The time consumed by the modified MNF transforma-als according to their histogram distribution, refine thetion was less than 22 seconds. The experimental resultsunmixed images by using a morphological operation andshowed that the computational efficiency of the modi-obtain the final endmembers.fied method was almost five times faster than that of thetraditional algorithm and that the computational per-4 Experiments and Resultsformance of the traditional method was greatly improvedby introducing the spectral subspace-processing4.1 Data Descriptionscheme.The data set used in this paper is an AVIRIS dataIn the next experiments, the ICA mixture modelset of a mixed agriculture/ forestry landscape acquiredwas used to resolve the mixing matrix and to obtain thein the Indian Pine Test Site ，Northwestern Indiana ,unmixing materials. After 69 steps of iteration with theJune 1992. The data set was composed of 220 spectralinfomax algorithm, the mixing matrix was obtained ，bands acquired in the 0.4 ~ 2.5 μum wavelength range.and the unmixing materials were obtained by inverseRemoving those bands that are corresponding to the wa-mapping with the mixing matrix. In order to reduce theter absorption regions, low SNR and bad bands, 200disadvantageous effect of the residual， each unmixingbands remain available. A scene of 145 x 145 pixels inmaterial was segmented according to its histogram dis-size was selected for our experiments, and it includestribution and the final endmembers shown in Fig. 4some ground cover， such as corn， soybean, wood,were obtained. To analyze the experimental results,wheat, grass, man-made materials, etc. The groundthe available ground truth map was referenced and usedtruth map is available in almost the whole scene. Theto identify which ground cover materials matched theground sampling distance of the hyperspectral images isendmembers obtained by the proposed algorithm. From25 m. Fig. 3 shows band 6 of the sensor.the result images shown in Fig. 4, some materials werebetter identified， such as woods， windrowed hay,grass/ trees，grass/ pasture, wheat, man- made objects ,etc.5 ConclusionsIn this paper， an unsupervised linear spectralmixture analysis algorithm was proposed for SMA in hy-perspectral images. This algorithm includes two keyFig.3 The hyperspectral image from band 6procedures: the modified MNF transformation and the4.2 Experimental ResultsICA mixture model. Through introducing the spectralIn the subspace partition of hyperspectral data,subspace-processing scheme, the modified MNF trans-the correlation threshold was set to be 0. 5. Thus fullformation is performed in data subspace. This process-data with 200 bands was partitioned into five sub-ing greatly improves the computational efficiency, andspaces: 1 ~35, 36~38 39 ~80, 81 ~ 103，104 ~the time consumed by the modified method is almost200. Fig. 1 describes the correlation between neighbor-one ffth of the traditional MNF transformation. Thening bands according to the correlation matrix of AVIRISthe I-unmix hyperspectraldata with 200 bands. From this figure, it was founddatathat the correlation between neighboring bands ap-withouMYHC N M H Ghistogram segmenta-peared to have blocking features.tion and morphological operation are used to obtain theln our experiments, the noise covariance estima-final endmembers, according to the histogram distribu-tion with shift difference method was performed on ation of the resulting images. Thus, the proposed algo-..475.Journal of Harbin Institute of Technology (Nev Series), Vol. 12, No. s, 2005rithm realizes the unsupervised linear spectral unmixingmote Sensing, 2000, 38(5): 2346 - 2360.of hyperspectral data with higher computational effi-[5] CHANG Cheng 1, ZHAO X L, ALTHOUSE MLC,etciency.al. Least squares subspace projection approach to mixedpixel casificatin for hyperspectral images [J]. IEEETransactions Geoscience and Remote Sensing, 1998, 36(3): 898 -912.[6] ZHANG Junping, ZHANG Ye, zou Bin, et al. Fusionlasification of hyperspectral image based on adaptivesubspace decompsition [ A ]. Intemational Conferenceon Image Processing ( ICIP2001) [C]. Thessaloniki:IEEE, 2001.[7GU Yanfeng, ZHANG Ye, ZHANG Junping. A kernelbased nonlinear subspace projection method for reduction(a) Soybeans(b) Comof bhyperspectral image dimensionality[ A]. IntermationalConference on Image Processing ( ICIP2002) [C]. Ro-chester, 'New York: lEEE, 2002.[8] BAYLISS J, GUALTIERI J A, CROMP R F. Analyzinghyperspectral data with independent component analysis[ EB/OL]. Proceedings of SPIE Applied Image and Pat-tem Recognition Workshop, 1997( htp://ww. cs. roch-ester. edu/ w bayliss spectral/spectral. html) .[9] TUTM, HUANG PS, CHEN P Y. Blind separation ofspectral signatures in hyperpectral imagery[J]. IEE(0) Hay-windrowed(d) GrasspasturesProe Vis. Image Signal Processing, 2001, 48(4): 217-226.[10] LEE Te- won, LEWICKI M s. Unsuperised imageclassification, segmentation, and enhancement usingICA mixture models [J]. IEEE Transactions on ImageProcessing, 2002, 11(3): 270-279.[11] MAKEIGS, BELL A J, JUNG T P, et al. Independentcomponent analysis of electroencephalographic data [J].Advances in NeuralInformnation Processing Systems,1996, 8: 145- 151.(e) Road, Stone-strel towers() Wheat[12] GREEN A A, BERMAN M, SWTrZER P,etal. AFig.4 The final endmembers after post-processingtransformnation for ordering mulispectral data in terms ofimage quality with implications for noise removal [J].IEEE Trans Greosei Renote Sensing, 1998. 26<1): 65References:[1] CHANGCI, DUQ, TU T M, et al. A joint band priori- .[13] GORDON C. A generalization of the maximum noisetization and bandccorelation approach to band selec-fraction transform [ J]。IEEE Trans Ceosci Remotetion for hyperspectral image lasfcation [ J]. IEEESensing, 2000, 38(1): 608 - 610.Trans Geosci Remote Sensing, 1999, 37(6): 2631 - 。[14] ROGER R E. A faster way to compute the noise adjus-_2640.ted principal components transform matrix [J]. IEEE[2SERPICOS B, BRUZZONE L A new search algorithmTrans Geosei Remote Sensing, 1994, 32(6): 1194 -for feature selection in hyperspectral remoteemote sensing ima-1 196.ges [J]. IEEE Trans Geosei Remote Sensing, 2001 , 39[15] PEARLMUTTER BA, PARRAL C. A context - sensi-(7): 1360 - 1367. .tive generalization of ICA [ EB/OL]. Proe. ICONIP'[3KFSHAVA N, MUSTARD J F. Spectral unmixing [J].96，1996. Available at_ fp://fp. cnl. salk. edu/IEEE Signal Prcessing Magpaine, 2002, 19(1): 44 -57.pub/bap/iconip - 96 - cica. ps. gz.[4] BROWN M, LEWIS H C, CUNNS R.Linear spectral[16] HYVARINEN A. Fast and robust fixed-pint algorithmmixture models and support vector machines for remotefor independent component analysis [ J]. IEEE TransNeural Network, 1999, 10(3): 626 -634.sensing [小]. IEEE Transactions on Geoscience and Re-中国煤化工MHCNM HG