GRAPH REGULARIZED COUPLED SPECTRAL UNMIXING FOR CHANGE DETECTION
Naoto Yokoya1 and Xiaoxiang Zhu2,31

Department of Advanced Interdisciplinary Studies, University of Tokyo, Japan
Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Germany
3
Chair of Remote Sensing Technology, Technische Universität München (TUM), Germany
2

ABSTRACT

Index Terms— Coupled spectral unmixing, graph regularization, change detection, multitemporal analysis
1. INTRODUCTION
Change detection has been an important task in multitemporal remote sensing data analysis. Spaceborne hyperspectral
imaging systems are promising for monitoring and characterising the earth’s environment on a global scale owing to its
advanced ability for accurate identification and classification.
Several spaceborne hyperspectral imaging missions, such as
the environmental mapping and analysis program (EnMAP)
and hyperspectral imager suite (HISUI), are going to be operational in the near future. An analysis of multitemporal hyperspectral data is expected to be useful for understanding
dynamic changes of the earth’s surface in details for a wide
range of applications, for example, disaster management, forest monitoring and precision agriculture. However, it raises
the demand for sophisticated algorithms to deal with the highdimensional nature of the data as well as multitemporal data…
Endmember
extracon0.7
Manmade 1
Water
Manmade 2

0.6Vegetation
Manmade 3
Soil

Coupled spectral unmixing

0.5

ReflectanceThis paper presents a methodology of coupled spectral unmixing for multitemporal hyperspectral data analysis. Coupled spectral unmixing simultaneously extracts the sets ofspectral signatures of endmembers and respective abundancemaps from multiple spectral images with differences in observation conditions and sensor characteristics. The problem
is formulated in the framework of coupled nonnegative matrix factorization. A graph regularization that reflects spectral
correlation between two images on abundance fractions isintroduced into the optimization of coupled spectral unmixing to consider temporal changes of the earth’s surface. An
alternating optimization algorithm is investigated using themethod of Lagrange multipliers to guarantee a stable convergence. The proposed method was applied to dual-temporal
Hyperion images taken over the Fukushima Daiichi nuclear
power plant. Experimental results showed that the proposed
method can extract essential information on the earth’s surface in a data-driven manner beyond multitemporal data
modality.

me ( )

0.4
0.3
0.20.10400600

800 1000 1200 1400 1600 1800 2000 2200 2400

Wavelength (nm)…

Fig. 1. Concept of coupled spectral unmixing.

modality due to the effect of atmospheric conditions, sensor
calibration and ground conditions, etc.
Supervised classification methods have been commonlyused for class-specific change detection using multitemporal spectral images [1]; however, they usually require a large
amount of training samples to obtain an accurate classifier,which often cost high in terms of human and economic resources. Spectral unmixing has been actively developed asa key technology in hyperspectral image processing to retrieve spectral signatures of endmembers and their fractional
abundance maps in a data-driven manner. When applied tomultitemporal data analysis, spectral unmixing is expected to
be useful not only for class-specific change detection [2] but
also for within-class change detection. However, few studies
have been carried out so far on multitemporal spectral unmixing that uses spectral correlation between multitemporal spectral images to improve the unmixing accuracy and robustness
against multitemporal data mismatch.
Multitemporal coupled spectral unmixing is defined as aprocess that simultaneously estimates the sets of spectral signatures of endmembers and abundance maps from multitemporal hyperspectral data in a data-driven manner. The benefits of multitemporal hyperspectral image analysis based on
spectral unmixing are threefold: 1) data-driven information
extraction of the earth’s surface owing to rich spatial-spectral-

temporal data; 2) class-specific change detection; 3) withinclass change detection. This paper presents a methodologyof multitemporal coupled spectral unmixing for land-coverchange detection. Coupled spectral unmixing for multitemporal hyperspectral data is formulated in the framework of coupled nonnegative matrix factorization [3] using graph regularization on spectral correlation between multitemporal data.
Experiments are conducted on the analysis of dual-temporal
Hyperion images taken over the Fukushima Daiichi nuclear
power plant in 2012 and 2014 to visualize the changes of landcover classes.
2. METHODOLOGY

The following update rules can be derived by the method ofLagrange multipliers

X = EAT + N(1)

where E ∈ RB×M is the spectral signature matrix, with M
being the number of endmembers, A ∈ RP ×M is the abundance matrix, N ∈ RB×P is the residual, and the operator ()Tdenotes the transposition operation. Researchers have studied
many models and algorithms based on geometrical, statistical,
and sparse regression-based approaches searching for robust,stable, tractable, and accurate unmixing [4].
Over the last decade, nonnegative matrix factorization(NMF) algorithms [5, 6] have emerged as useful unmixing
methods because they overcome the difficulty of absence ofpure pixels with straightforward implementation and mathematical flexibility to consider several constraints. Lu etal. presented a graph regularized L1/2 -NMF (GLNMF) thattakes into account sparse characteristic and the intrinsic manifold structure of the data [7]. With a matrix W ∈ RP ×P
representing the weight matrix of the graph, the optimization
is given bymin ∥X − EAT ∥2F + α
E,A+βP∑P ∑M∑1/2Ajm

j=1 m=1∥aj − al ∥22 Wjl(2)s.t. E ⪰ 0, A ⪰ 0,M∑

Ajm = 1m=1where aj = [aj1 , ... , ajM ]T , ∥ · ∥F denotes the Frobenius
norm, and the symbol ⪰ denotes inequality component-wise.

(4)

D ∈ RP ×P denotes a diagonal matrix whose entries are column sums of W. There are many choices to define the weightmatrix W on the graph, for example, 0-1, heat kernel, anddot-product weighting. When both spatial and spectral distances are considered to construct the graph, each element of
W is given by(Wjl = e−

∥Xj −Xl ∥2d(j,l)
2σd −σs)(5)where d(j, l) denotes the Euclidean distance between the j thand its neighboring lth pixel, and σd and σs are the parameters
for spatial and spectra distances, respectively.

2.2. Graph-regularized coupled sparse nonnegative matrix factorization
We extend GLNMF to multitemporal spectral unmixing by
integrating the graph regularization on the temporal changesof spectra into two unmixing problems in a coupled basis.
We assume that when the spectrum has changed, the abundance fractions have also changed, and vice versa. In thispaper, we consider coupled spectral unmixing of two spectral images denoted by X1 ∈ RB1 ×P and X2 ∈ RB2 ×P .
The two images are geometrically co-registered with the sameground sampling distance but the band numbers can be different to use multiple sensors. Consider a weight matrix Wt ∈RP ×P with each element Wtjl corresponds to the similarity between X1j and X2l . Graph-regularized coupled sparse
nonnegative matrix factorization (GCSNMF) can be formulated as2
T 2
minE1 ,A1 ,E2 ,A2 ∥X1 −E1 AT
1 ∥F +∥X2 −E2 A2 ∥F)∑P ∑M (1/21/2+α j=1 m=1 A1 jm +A2 jm+β∑P+γ∑Pj,l=1(∥a1 j −a1 l ∥22 W1 jl +∥a2 j −a2 l ∥22 W2 jl )22j,l=1 ∥a1 j −a2l ∥2 Wt jl +δ∥RE1 −E2 ∥F

s.t. E1 ⪰0, A1 ⪰0, E2 ⪰0, A2 ⪰0,

j,l=1

(3)

(XT E+βWA)jmAjm ← Ajm AET E+ α A1/2 +βDA()jm
2

2.1. Spectral unmixing based on nonnegative matrix factorizationA linear spectral mixture model is commonly used for unmixing problems owing to its physical effectiveness and mathematical simplicity. A spectrum at each pixel is assumed to bea linear combination of several endmember spectra. Therefore, a hyperspectral image X ∈ RB×P , with B bands and Ppixels, is formulated as(XA)im
Eim (EAT A)
im

Eim ←∑m A1 jm =1,∑(6)

m A2 jm =1

where E1 ∈ RB1 ×M1 and A1 ∈ RP ×M1 are the signatureand abundance matrices of the first data X1 , E2 ∈ RB2 ×M2and A2 ∈ RP ×M2 are those of the second data X2 , andR ∈ RB2 ×B1 is the relative spectral response function. In addition to the sparsity constraint and the spatial-spectral graphregularization for each unmixing, the term with the parameterδ is the temporal graph regularization. The update rules canbe obtained by the method of Lagrange multipliersT

1 A1 +δR E2 )im
E1im (E1(XAT A1 +δRT RE1 )im

E1im ←A1jm ←(7)

1A1jm ((XT1 E1 +βW1 A1 +γWt A2 )jm
1/2αA1 ET
1 E1 + 2 A1)+βD1 A1 +γDt1 A1E2km ←(X2 A2 +δRE1 )km
E2km (ET2 A A2 +δE2 )km

A2jm ←TXT
2 E2 +βW2 A2 +γWt A1 jm)
2jm (1/2α+βD2 A2 +γDt2 A2A2 ET
2 E2 + 2 A2(8)
jm

(9)

2A

()

(10)
jm

where Dt1 and Dt2 denote diagonal matrices whose entries
are column and row sums of Wt , respectively. Wt is designed as a diagonal matrix using the spectral angle distancefor the spectral similarity measurement. The abundance sumto-one constraint is implemented using a method given in [9].
After the initialization of all the matrices, they are alternately optimized by Equations (7)–(10). In this paper, wepresent a specific implementation of GCSNMF. Firstly, twoimages are obtained by the same sensor, and thus δ is set to 0.
Secondly, the two images mostly include the same endmembers with spectral mismatch. Thirdly, the data is assumed tocontain all pure pixels of analyzed materials.
The data matrix X1 is considered as the reference data andthe endmember matrix E1 is estimated by vertex componentanalysis (VCA) [8]. The endmember matrix E2 is initializedas [E1 r], where r ∈ RB is a mean vector of the residual between the two images, i.e., X2 −X1 , with all negative components replaced by a very small positive value. The abundance
matrices A1 and A2 are initialized by the fully constrained
least-squares (FCLS) method [9]. During the optimization ofGCSNMF, A1 and A2 are optimized by Equations (8) (10),respectively. For the spectral signatures, only the last column
vector of E2 is updated as a virtual endmember signature thatexplains spectral mismatch between the two images. Classspecific changes of land cover can be visualized by comparing multitemporal abundance maps. All the parameters in thecost function are empirically set taking a balance of the regularization terms.
3. EXPERIMENTThe proposed method is applied to dual-temporal Hyperion
images taken over the Fukushima Daiichi nuclear power planton 29 April 2012 and 25 May 2014. The hyperspectral data
cubes with composite RGB images are shown in Fig. 2. Atsensor reflectance data that include atmospheric effects areused with 197 effective bands. The first image taken on 29
April 2012 is used for the reference to extract endmembers.
In this experiment, totally six endmembers are extracted by
VCA and labeled by human as three types of man-made materials, vegetation, soil, and water. When applied to at-surface
reflectance data, a spectral dictionary is useful to make thislabeling process automatic. The proposed coupled spectral

1 km

Fukushima Daiichi
nuclear power plantFig. 2. Hyperion images taken over the Fukushima Daiichi
nuclear power plant on (left) 29 April 2012 and (right) 25
May 2014, respectively.

unmixing method is compared with the conventional individual spectral unmixing methods, i.e., FCLS and GLNMF. Notethat the same endmember signatures are used for all methods.
Fig. 3(a) shows the unmixing results on the first image obtained by the three methods, whereas Fig. 3(b) shows those on
the second image. The last row images are composite RGB
images using abundance maps of man-made materials, vegetation, and soil for red, green, and blue, respectively. Fig.
3(c) shows the difference maps of abundances with RGB images visualizing class-specific changes. For example, magenta, green, orange, and purple indicate the changes fromvegetation to man-made, from soil to vegetation, from soil toman-made, and from vegetation to soil, respectively. In theFCLS results, large abundance errors of water appear in theland area as a result of the strict sum-to-one constraint. In theresults obtained by FCLS and GLNMF, there are too muchabundance increase in vegetation caused by data mismatch
due to atmospheric effects. The proposed method mitigatedthis error owing to the temporal constraint. In addition, theresults obtained by the individual spectral unmixing methods
show too much abundance increase of man-made materials in
the Fukushima Daiichi nuclear power plant and the urban area
located on the northwest part of the power plant, although
they have just partly changed rather than largely. The proposed method shows stable and reasonable results on theseareas. Accordingly, the proposed method seems to achieve
the accurate spectral unmixing and its robustness against temporal data mismatch. In the RGB image of abundance difference obtained by GCSNMF, trees that have been cut down
to increase the number of water storage tanks in and aroundthe power plant can be clearly visualized as magenta and purple pixels. The increase of vegetation may be because of thegrowth of weeds owing to seasonal difference as well as theyearly expansion of the abandoned farmland in the restrictedarea.

FCLS

GLNMF

GCSNMF

FCLS

GLNMF

(a)

(b)

GCSNMF

FCLS

GLNMF

GCSNMF(c)Fig. 3. Unmixing results of data acquired on (a) 29 April 2012 and (b) 25 May 2014, and (c) their difference images obtainedby FCLS, GLNMMF, and GCSNMF from left to right. Abundance maps of man-made materials, vegetation, soil, and water,
and composite RGB images using abundance maps of man-made materials, vegetation, and soil for red, green, and blue,
respectively, are shown from top to bottom. The figure can be best viewed with 200% zoom.
4. CONCLUSION
This paper present a coupled spectral unmixing method for
multitemporal hyperspectral data analysis. The proposed
method introduces a graph regularization on spectral changes
into the optimization of multitemporal spectral unmixing ina coupled basis in order to improve the unmixing accuracy
and robustness against multitemporal data mismatch. The
proposed method has been applied to dual-temporal Hyperion images taken over the Fukushima Daiichi nuclear power
plant and successfully visualized class-specific changes ofland cover. Coupled spectral unmixing is expected to be apowerful tool for multitemporal spectral data analysis enabling data-driven extraction of essential information on theearth’s surface that cannot be obtained from a single image. Our future work includes investigation on numericalvalidation and application to multisensor data.
5. REFERENCES[1] G. Camps-Valls, L. Gmez-Chova, J. Muoz-Mar, J. L. Rojo-lvarez, and M.
Martnez-Ramn, “Kernel-based framework for multitemporal and multisource
remote sensing data classification and change detection,” IEEE Trans. Geosci.
Remote Sens., vol. 46, no. 6, pp. 1822-1835, 2008.

[2] D. Lu, M. Batistella, and E. Moran, “Multitemporal spectral mixture analysis
for Amazonian land-cover change detection,” Can. J. Remote Sens., vol. 30, no.
1, pp. 87–100, 2004.
[3] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix factorization
unmixing for hyperspectral and multispectral data fusion,” IEEE Trans. Geosci.
Remote Sens., vol. 50, no. 2, pp. 528–537, 2012.
[4] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, andJ. Chanussot, “Hyperspectral unmixing overview: geometrical, statistical, andsparse regression-based approaches,” IEEE J. Sel. Topics Appl. Earth Observ.,
vol. 5, no. 2, pp. 354–379, 2012.
[5] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix
factorization,” Nature, vol. 401, pp. 788–791, 1999.
[6] D. Cai, X. He, J. Han, and T. S. Huang, “Graph regularized non-negative matrix
factorization for data representation,” IEEE Trans. Pattern Anal. Mach. Intell.,vol. 33, no. 8, pp. 1548–1560, 2010.
[7] X. Lu, H. Wu, Y. Yuan, P. Yan, and X. Li, “Manifold regularized sparse NMF
for hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 51, no. 5,pp. 2815-2826, 2013.
[8] J. M. P. Nascimento and J. M. B. Dias, “Vertex component analysis: A fastalgorithm to unmix hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol.
43, no. 4, pp. 898–910, 2005.
[9] D. C. Heinz and C.-I. Chang, “Fully constrained least squares linear spectral
mixture analysis method for material quantification in hyperspectral imagery,”IEEE Trans. Geosci. Remote Sens., vol. 39, no. 3, pp. 529–545, 2001.