GRAPH-REGULARIZED COUPLED SPECTRAL UNMIXINGFOR MULTISENSOR TIME-SERIES ANALYSIS
Naoto Yokoya1,2,3 , Xiao Xiang Zhu1,2 , and Antonio Plaza4
1

2

Remote Sensing Technology Institute (IMF), German Aerospace Center (DLR), Germany
Signal Processing in Earth Observation (SiPEO), Technische Universität München (TUM), Germany
3
Department of Advanced Interdisciplinary Studies, University of Tokyo, Japan
4
Department of Technology of Computers and Communications, University of Extremadura, Spain
ABSTRACTA new methodology that solves unmixing problems involving
a set of multisensor time-series spectral images is proposedin order to understand dynamic changes of the surface at a
subpixel scale. The proposed methodology couples multiple
unmixing problems via regularization on graphs between themultisensor time-series data to obtain robust and stable unmixing solutions beyond data modalities owing to different
sensor characteristics and the effects of non-optimal atmospheric correction. A synthetic dataset that includes seasonal
and trend changes on the surface and the residuals of nonoptimal atmospheric correction is used for numerical validation. Experimental results demonstrate the effectiveness ofthe proposed methodology.
Index Terms— Coupled spectral unmixing, multisensor
data fusion, time-series analysis, change detection
1. INTRODUCTION
A synergetic use of hyperspectral and multispectral images
is important for upcoming spaceborne imaging spectroscopy
missions [1]. Spaceborne imaging spectroscopy enables monitoring dynamic processes of the surface in detail on a global
scale; however, revisit time is limited compared to multispectral imaging satellites. Understanding of the dynamics on thesurface can be improved by synergistically analyzing a timeseries of spaceborne hyperspectral and multispectral images.
In recent years, particular attention has been paid to spectral unmixing of multitemporal hyperspectral data owing to its
ability to detect class-specific changes at a subpixel level [2]–[11]. Previous research has shown the effectiveness of multitemporal spectral unmixing in a wide range of applications[2]–[6]. Current research is focused on applying spectrallibrary-based unmixing methods to multitemporal unmixing
This research has been partly supported by Alexander von Humboldt
Fellowship for postdoctoral researchers, Helmholtz Young Investigators
Group “SiPEO” (VH-NG-1018, www.sipeo.bgu.tum.de), and Japan Society
for the Promotion of Science (JSPS) KAKENHI 15K20955.

problems [5, 9]. A data-driven approach has also been developed in [10, 11].
To the authors’ knowledge, very few publications can befound in the literature that investigate time-series spectral unmixing from the viewpoint of synergistically exploiting hyperspectral and multispectral images. Hyperspectral and multispectral image fusion has received great attention to generate high-resolution hyperspectral data [12]; however, studies
on time-series analysis are still lacking. To tackle spectral
unmixing of a multisensor time series, it is challenging to determine how to use hyperspectral and multispectral images incomplementary fashion. Furthermore, most of the previous
studies on multitemporal hyperspectral unmixing do not take
into account possible data mismatches due to non-optimal atmospheric correction, even though it always appears in timeseries spectral data [13]–[15].
In this paper, a novel methodology for spectral unmixingof multisensor time-series data is proposed. The objective isto understand dynamic changes on the surface at a subpixel
scale from multisensor time-series data, beyond data modalities, due to different sensor characteristics and the effects of
non-optimal atmospheric correction. Each spectral unmixing
problem is solved cooperatively with those related to otherimages that have similar surface conditions to obtain a robust,
stable, and accurate solution. Multiple unmixing problems
are coupled via regularization on graphs between neighborhood images.
The remainder of the paper is organized as follows: Section II introduces our methodology. Experimental results on
synthetic data are presented in Sections III. Section IV concludes the paper.
2. METHODOLOGY
To obtain robust and stable unmixing beyond data modalities
due to different spectral characteristics and non-optimal atmospheric correction, we introduce coupled spectral unmixing
for time-series analysis of multisensor spectral data, whichsolves an unmixing problem of a single image cooperatively

of coupled spectral unmixing can be formulated as

K+1

1 X
kYk − Ak XTk k2F
X,Z 2

min

k=1

K

K+1

P

1 X X X+ α((G)pq (G)qp kXjp − Xlq k22
2 p=1 q=p+1

(1)

j,l=1

Hyperspectral data
Multispectral data+((G)qp − (G)pq )kXjp − Zlq k22 )(Wpq )jl )s.t. Xk  0, Xk = Zk ,Fig. 1. Manifold of multisensor time-series spectral datawith those of its neighboring images (in the time-series hierarchy) by adopting an intrinsic manifold of the time-series
dataset. Fig. 1 illustrates the manifold of multisensor timeseries spectral data. Each point represents a single image.
Owing to changes on the surface, the multisensor time-series
dataset forms the manifold, where neighborhood images have
similar conditions of the surface. Coupled spectral unmixingof multisensor time-series data is formulated as a decentralized optimization problem for multi-agent systems [16],where each spectral unmixing is solved via a local optimization problem based only on information concerning otherspectral unmixing problems in its neighborhood.
The multiple spectral unmixing problems are coupled via
regularization on local graphs. In a local graph, a target image and its neighbors are regarded as nodes and connected byedges. Local graphs can be mixtures between undirected anddirected graphs, depending on the difference in the spectral
resolution between the two images. Undirected graphs areused between images with similar spectral resolution so that
local spectral unmixing problems mutually influence eachother. Hyperspectral imaging enables more accurate spectral
unmixing than multispectral imaging. If the target image is amultispectral data set, hyperspectral neighboring images canimprove the accuracy of the spectral unmixing for the target
image; however, the reverse is not necessarily the case. Therefore, directed graphs, where the edges have a direction fromhigher-spectral-resolution data to lower-spectral-resolution
data in the manifold, are also considered (Fig. 1).
Let us consider K + 1 images {Yk ∈ RBk ×P }K+1
k=1 , i.e.,a target image and its K-nearest neighbors. The numbers ofspectral bands satisfy B1 ≤ B2 ≤ ... ≤ BK+1 without lossof generality. For each spectral unmixing, L1/2 sparsity regularization on abundances is adopted []. It is supposed thatthe spectral signature matrix is obtained either from a spectral library or from the data itself. The optimization problem

where X = [XT1 , ..., XTK+1 ]T , Xk ∈ RP ×M denotes theabundance matrix of the kth data (Yk ). The penalty termwith the parameter α is the graph regularization, which places
the restriction such that, if two spectral signatures are similar,
the abundance vectors are also similar to each other. Wpq ∈RPp ×Pq is the weight matrix on the graph between the pth
and qth data, and (Wpq )jl represents the closeness of two
points: the jth column of Yp and the lth column of Yq . Wpq
is designated by the spectral similarity measurement using a
distance metric, such as the Euclidean distance or spectral angle distance (SAD). G ∈ RK+1×K+1 is the adjacency matrix that represents the local graph between K + 1 images.
Gpq = Gqp = 1 if the pth and qth images are connected by
the undirected edge, whereas Gpq = 0 and Gqp = 1 if they
are connected by the directed edge from the qth image to thepth image. Z = [ZT1 , ..., ZTK+1 ]T is auxiliary variables that
implement the directed-graph regularization.
The formula in (1) can be solved by alternating optimization after the initialization of X and Z using a conventional
method for individual spectral unmixing, such as the constrained least squares (CLS) method [17]. The alternating
update rules (2)–(4) are obtained by the augmented Lagrange
method. Θk  0 denotes the augmented Lagrangian parameters. Dpq 1 ∈ RPq ×Pq and Dpq 2 ∈ RPp ×Pp denote diagonal
matrices whose entries are column and row sums of Wpq ,respectively. The simplest way to construct Wpq is to calculate the spectral similarity between only a pair of spatiallycorresponding pixels. In this case, Wpq = Dpq1 = Dpq2 if
Pp = Pq is satisfied. We use this technique to save computational costs. The average of non-zero elements in W is used
to find the K-nearest neighbors and construct G, which is referred to as manifold-based coupling. An alternative means ofconstructing G is the use of sequential coupling that uses temporally neighboring images. Note that the problem (1) does
not include the abundance sum-to-one constraint. To satisfy
the abundance sum-to-one constraint, a method given in [17]is adopted.
The proposed algorithm is summarized in Algorithm 1.

(Xk )jm ← (Xk )jm

Pk−1
PK+1
T
p=1 Gpk Gkp Wpk Xp +q=k+1 ((G)kq (G)qk Wkq Xq +((G)qk −(G)kq )Wkq Zq ))+Θk )jm
Pk−1
PT(Xk Ak Ak +α( p=1 (G)pk (G)kp Dpk2 Xk + K+1
q=k+1 ((G)kq (G)qk +(G)qk −(G)kq )Dkq1 Xk ))jmT
(Yk
Ak +α(Tα( p=1 ((G)kp −(G)pk )Wpk
Xp +ρXk )jmjm (α Pk−1 ((G) −(G) )DZ+Θk +ρZk )jm
kp
pk
pk2
k
p=1

(Zk )jm ← (Zk )

(2)

Pk−1

(Θk )jm ← (Θk )jm + ρ(Zk − Xk )jm

Algorithm 1 Graph-regularized coupled spectral unmixing
for multisensor time-series data
Input: Y, A
Output: X
1: for i = 1 to N do
2:
calculate W between all combinations between Yi
and the others
3:
find K nearest neighbors of Yi
4:prepare G
5:initialize {Xk }k∈Vi by individual spectral unmixing(Vi is a set of the ith target and K neighbor images)
6:
solve (1) by alternating updates of (2)–(4)7:update {Xi }8: end for
9: return {Xi }N
i=1

3. EXPERIMENT
3.1. Synthetic multisensor time-series spectral images
A set of synthetic multisensor time-series spectral images wasused to numerically validate the proposed methodology. The
synthetic dataset was simulated using the abundance maps of
the Fractal 1 synthetic data [18]. The Fractal 1 image is composed of nine endmembers; their abundance maps were generated from fractal patterns with the size of 100 × 100 pixels.
To simulate realistic changes on the surface, seasonal changes
and trend changes were considered. Nine endmembers wereselected from the U.S. Geological Survey (USGS) spectral
library: grass, dry grass, oak, soil, melting snow, water, asphalt, green house, and concrete. Seasonal changes were simulated for grass, dry grass, oak, soil, and melting snow using
periodic functions, and trend changes are simulated for green
house, and concrete using a monotonic increasing function.
The spectral response functions of AVIRIS (bands 1–32,
36–96, 100–160, and 163–224) and Landsat-8 (bands 1–8)
were used to simulate hyperspectral and multispectral images,
respectively. Additive Gaussian noise with a signal-to-noise
ratio (SNR) of 100 was added to each temporal image. Nonoptimal atmospheric correction was simulated using a Gaussian random vector for the residual gains. The revisit cycles of
the hyperspectral and multispectral imagers were set to 16 and27 days, respectively, which correspond to those of Landsat-8
and EnMAP. Two scenarios are considered for data acquisi-

(3)
(4)

tion: 1) full acquisition, that assumes a clear sky for all images, and 2) realistic acquisition, that assumes a clear sky for
30% of the observations. The total observation period is set tofive years and therefore the total number of images is 68 and
115 for hyperspectral and multispectral images, respectively.
Composite color images of five examples are shown in Fig.
2(a). The corresponding reference abundances of grass andconcrete are presented in the first row images of Fig. 2(b)(c).
3.2. Experimental results
We use three methods: 1) individual spectral unmixing (ISU)based on the CLS method; 2) coupled spectral unmixing withsequential coupling (CSU 1); and 3) coupled spectral unmixing with manifold-based coupling (CSU 2). To investigate
the efficacy of using multiple sensors, these methods are performed on three different datasets: 1) only hyperspectral images; 2) only multispectral images; and 3) hyperspectral andmultispectral images. A is given, K = 2, and α = 1.
Table 1 shows the average root-mean-squareerrors (RMPN
SEs) of abundances, defined by N1 i=1 kXi − X̃i kF . For
the realistic data acquisition scenario, the average results of
10 Monte Carlo trials obtained after selecting 30% of theobservations are shown. Note that only the unmixing accuracy of multispectral images is shown for the third dataset,because only multispectral images can gain benefits, whereas
the results of hyperspectral images are the same as for the
first dataset. The coupled spectral unmixing algorithms outperform individual spectral unmixing for the first and third
datasets. The unmixing accuracy of multispectral data ishighly improved owing to the support of neighboring hyperspectral images. Coupled spectral unmixing based on
manifold-based coupling outperforms that based on sequential coupling in many cases, particularly in the realistic dataacquisition scenarios. This suggests that manifold-based
coupling leads to more accurate and robust unmixing results than sequential coupling. Coupled spectral unmixing
with sequential coupling performs well in ideal scenarios
because sequential neighbors with high temporal resolution
represent neighbors in the manifold of time-series data well.
In contrast, if data acquisition is limited to realistic scenarios, manifold-based coupling is required because sequential
neighbors do not always represent neighbors in the manifold.
Figs. 2(b)(c) present reference abundance maps of grassand concrete for the five examples of the third dataset with

Table 1. Average RMSEs of abundances obtained before atmospheric normalization from three datasets of synthetic time-series
spectral images, i.e. only hyperspectral (HS) data (68 images), only multispectral (MS) data (115 images), and both HS andMS data (183 images) using ISU, CSU 1, and CSU 2. RMSEs for only MS data are shown for the last dataset.
DatasetAcquisition
ISU
CSU 1
CSU 2

HS
Ideal
0.028181
0.025268
0.025552

Realistic
0.028526
0.026195
0.025896

MS
Ideal
Realistic
0.11484 0.11719
0.11610.121100.11538 0.11840

MS (coupled with HS)
Ideal
Realistic
0.11484
0.117190.086309 0.092507
0.049207 0.052016

Reference

(c)

CSU 1

CSU 1

ISU

ISU

(b)

Reference

(a)

1

CSU 2

CSU 2

1

0

0

Fig. 2. (a) Color composites of ten images in the synthetic dataset and their corresponding abundance maps of (b) grass and(c) concrete. Abundances of reference and estimated abundances obtained by ISU, CSU 1, and CSU 2 are shown from top tobottom in (b) and (c).
one trial of the realistic acquisition scenario and the corresponding estimated abundance maps obtained by the threemethods. Serial day numbers of the subset are 433, 1009,
1313, 1473, and 1825, which are all simulated as observed by
the multispectral imager. In Fig. 2(b), the coupled spectral
unmixing methods shows better resemblance to the referencecompared to individual spectral unmixing. In particular, CSU
2 exhibits the best results. On the other hand, in Fig. 2(c),CSU 1 shows the best resemblance, whereas both CSU 1 andCSU 2 outperform individual spectral unmixing. These results imply that manifold-based coupling better captures seasonal changes and sequential coupling is more suitable fordetection of trend changes.

data modalities due to the different spectral characteristics
of imagers and the effects of non-optimal atmospheric correction. A single unmixing problem is cooperatively solved
with other unmixing problems in its neighborhood in the manifold of the time-series dataset. Multiple unmixing problems
are coupled via regularization on graphs between the local
images. Numerical validation was performed with synthetic
time-series data composed of hyperspectral and multispectral
images, including trend and seasonal changes on the surface
and residual gains simulating non-optimal atmospheric correction. The experiment demonstrated the effectiveness of thecoupled spectral unmixing algorithm compared to individual
spectral unmixing. A generalization of the proposed methodology and its application to real data is noteworthy [19].

4. CONCLUSIONS
5. REFERENCES
A coupled spectral unmixing methodology was proposed for
multisensor time-series analysis. The proposed methodology provides robust and stable unmixing solutions beyond[1] L. Guanter, H. Kaufmann, K. Segl, S. Foerster, C. Rogass, S. Chabrillat, T. Kuester, A. Hollstein, G. Ross-ner, C. Chlebek, C. Straif, S. Fischer, S. Schrader, T.
Storch, U. Heiden, A. Müller, M. Bachmann, H. Mühle,
R. Müller, M. Habermeyer, A. Ohndorf, J. Hill, H. Buddenbaum, P. Hostert, S. van der Linden, P. J. Leitão, A.
Rabe, R. Doerffer, H. Krasemann, H. Xi, W. Mauser, T.
Hank, M. Locherer, M. Rast, K. Staenz, and B. Sang,“The EnMAP spaceborne imaging spectroscopy mission for Earth observation,” Remote Sensing, vol. 7, pp.
8830–8857, 2015.
[2] D. Lu, M. Batistella, and E. Moran, “Multitemporal spectral mixture analysis for Amazonian land-coverchange detection,” Can. J. Remote Sens., vol. 30, no. 1,
pp. 87–100, 2004.
[3] D. B. Lobell and G. P. Asner, “Cropland distributions
from temp oral unmixing of MODIS data,” Remote Sens.
Environ., vol. 93, no. 3, pp. 412–422, Nov. 2004.
[4] M. A. Goenaga, M. C. Torres-Madronero, M. VelezReyes, S. J. Van Bloem, and J. D. Chinea, “Unmixing
analysis of a time series of Hyperion images over the
Guánica dry forest in Puerto Rico,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., vol. 6, no. 2, pp.
329–338, Apr. 2013.
[5] B. Somers and G. P. Asner, “Multi-temporal hyperspectral mixture analysis and feature selection for invasive
species mapping in rainforests,” Remote Sens. Environ.,
vol. 136, pp. 14–27, Sep. 2013.
[6] P. Du, S. Liu, P. Liu, K. Tan, and L. Cheng, “Subpixel change detection for urban land-cover analysis via
multi-temporal remote sensing images,” Geo-spatial Inf.
Sci., vol. 17, no. 1, pp. 26–38, 2014.
[7] Q. Du, L. Wasson, and R. King, “Unsupervised linear unmixing for change detection in multitemporal airborne hyperspectral imagery,” in Proc. IEEE Int. Workshop MultiTemp, 2005, pp. 136–140.
[8] A. Ertürk and A. Plaza, “Informative change detection
by unmixing for hyperspectral images,” IEEE Geosci.
Remote Sens. Lett., vol. 12, no. 6, pp. 1252–1256, Jun.
2015.
[9] A. Ertürk, M.-D. Iordache, and A. Plaza, “Sparse unmixing based change detection for multi-temporal hyperspectral images,” IEEE J. Sel. Topics Appl. Earth Observ. Remote Sens., in press.
[10] N. Yokoya and X. X. Zhu, “Graph regularized coupledspectral unmixing for change detection,” Proc. WHISPERS, Tokyo, Japan, June 2-5, 2015.
[11] S. Liu, L. Bruzzone, F. Bovolo, and P. Du, “Unsupervised multitemporal spectral unmixing for detecting
multiple changes in hyperspectral images,” IEEE Trans.
Geosci. Remote Sens., in press[12] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix factorization unmixing for hyperspectral andmultispectral data fusion,” IEEE Trans. Geosci. Remote
Sens., vol. 50, no. 2, pp. 528-537, 2012.
[13] A. Singh, “Digital change detection techniques using
remotely-sensed data,” Int. J. Remote Sens., vol. 10, no.
6, pp. 989–1003, 1989.
[14] C. Songh, C. E. Woodcock, K. C. Seto, M. P. Lenney,and S. A. Macomber, “Classification and change detection using Landsat TM data: when and how to correct
atmospheric effects?” Remote Sens. Environ., vol. 75,
pp. 230–244, 2001.
[15] P. Coppin, I. Jonckheere, K. Nackaerts, B. Muys, and E.
Lambin, “Digital change detection methods in ecosystem monitoring: a review,” Int. J. Remote Sens., vol. 25,
no. 9, pp. 1565–1596.
[16] S. L. Waslander, G. Inalhan and C. J. Tomlin, Decentralized Optimization via Nash Bargaining, 2004, Kluwer
Academic Press.
[17] D. C. Heinz and C.-I. Chang, “Fully constrained leastsquares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE
Trans. Geosci. Remote Sens., vol. 39, no. 3, pp. 529–
545, 2001.
[18] J. Plaza, E. M. T. Hendrix, I Garcı́a, G. Martı́n, A. Plaza,“On endmember identification in hyperspectral images
without pure pixels: A comparison of algorithms,” J.
Math. Imaging Vis., vol. 42, pp. 163–175, 2012.
[19] N. Yokoya, X. X. Zhu, and A. Plaza, “Multisensor coupled spectral unmixing for time-series analysis,” IEEE
Trans. Geosci. Remote Sens., submitted, 2016.