1430

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 52, NO. 2, FEBRUARY 2014

Nonlinear Unmixing of Hyperspectral Data Using
Semi-Nonnegative Matrix Factorization
Naoto Yokoya, Student Member, IEEE, Jocelyn Chanussot, Fellow, IEEE, and Akira Iwasaki

Abstract— Nonlinear spectral mixture models have recently
received particular attention in hyperspectral image processing.
In this paper, we present a novel optimization method of nonlinear unmixing based on a generalized bilinear model (GBM),
which considers the second–order scattering of photons in a
spectral mixture model. Semi-nonnegative matrix factorization(semi-NMF) is used for the optimization to process a whole
image in matrix form. When endmember spectra are given, theoptimization of abundance and interaction abundance fractions
converge to a local optimum by alternating update rules withsimple implementation. The proposed method is evaluated usingsynthetic datasets considering its robustness for the accuracyof endmember extraction and spectral complexity, and shows
smaller errors in abundance fractions rather than conventional
methods. GBM-based unmixing using semi-NMF is applied tothe analysis of an airborne hyperspectral image taken overan agricultural field with many endmembers, and it visualizes
the impact of a nonlinear interaction on abundance maps atreasonable computational cost.
Index Terms— Generalized bilinear model (GBM), nonlinear
unmixing, semi-nonnegative matrix factorization.

I. I NTRODUCTIONSPECTRAL unmixing is an important task for hyperspectral image interpretation. Many researchers have workedon this problem using a linear mixture model (LMM), in which
it is assumed that an observed spectrum is a linear combinationof several endmember spectra. The LMM is a simplified
spectral mixture model that considers only first-order scatteredphotons by neglecting multiple photon interactions. Although
LMM-based unmixing methods can retrieve physically meaningful results, nonlinearity in a spectral mixture model has
been pointed out in many works [1]–[12]. In recent years,
nonlinear unmixing for hyperspectral images has been receiving particular attention in remote sensing image exploitation.
Nonlinear spectral mixing occurs owing to multiple reflections
and transmissions from a surface [2]. Bioucas–Dias et al. [13]classified nonlinear unmixing into two scenarios: an intimate

Manuscript received October 16, 2012; revised February 21, 2013; accepted
February 22, 2013. Date of publication May 7, 2013; date of current version
December 12, 2013. This work was carried out in part under Contract with
the National Institute of Advanced Industrial Science and Technology andwas supported by the Hypanema Project Funded by Agence Nationale de laRecherche (ANR) under Grant ANR-12-BS03-0003.
N. Yokoya and A. Iwasaki are with the Department of Aeronautics
and Astronautics, University of Tokyo, Tokyo 153-0041, Japan (e-mail:yokoya@sal.rcast.u-tokyo.ac.jp; aiwasaki@sal.rcast.u-tokyo.ac.jp).
J. Chanussot is with the Grenoble Institute of Technology, Saint Martind’Hères cedex 38402, France (e-mail: jocelyn.chanussot@gipsa-lab.grenobleinp.fr).
Color versions of one or more of the figures in this paper are available
online at http://ieeexplore.ieee.org.
Digital Object Identifier 10.1109/TGRS.2013.2251349

mixture and a multilayered scene. The intimate mixture modelis based on the photometric model of Hapke [14], which
considers multiple scattering between different materials at
the particle level [6]. Close et al. [7], [8] applied Hapke’s
average albedo model to solve fully unsupervised nonlinear
unmixing in the case of intimate mixtures. In the multilayered scene, there are multiple interactions among scatters at
different layers, which often happen between vegetation andsoil [1]–[5]. The bilinear mixture model (BMM) considers thesecond-order scattering of photons between two distinct materials and introduces additional virtual endmembers formedby their product terms into the LMM [3]. The BMM is acommon model of nonlinear unmixing in multilayered scenes,
and its optimization has been studied with different constraints
by several groups [5], [9]–[11]. Chen et al. [16] analyzed
the impact of the collinearity effect [15] of virtual endmembers on bilinear spectral unmixing. When there is a highcorrelation between endmembers and second-order scattering
virtual endmembers and when hyperspectral data containsevere noise, the optimization based on [5] resulted in worse
abundance estimation than linear spectral mixture analysis.
Halimi et al. [11] introduced the generalized bilinear model(GBM) as an effective means of dealing with the underlying
assumptions in the BMM. The GBM method was applied tothe analysis of small images of synthetic and real hyperspectral
data with three endmembers and showed good results [11],[12]. When applied to the analysis of larger images in an
unsupervised manner with more endmembers, the optimization
becomes more challenging owing to the collinearity and thevarious possibilities of local minima.
Nonnegative matrix factorization (NMF) [17], [18], which
factorizes a nonnegative matrix into two nonnegative matrices,
has recently emerged as a useful method of solving theLMM-based unmixing problem [19]–[22]. Ding et al. [23]proposed a new variation in the theme of NMF, i.e., a seminonnegative matrix factorization (semi-NMF) that factorizes anonrestricted matrix into a nonrestricted matrix and a nonnegative matrix. In this paper, we present semi-NMF as a novel
optimization method for GBM-based unmixing. We evaluate
it using synthetic datasets and apply it to the analysis of a real
hyperspectral image with 10 endmembers. We also investigate
the effect of endmember extraction on nonlinear unmixing
and visualize the impact of second-order scattering effects onabundance maps.
The rest of this paper is organized as follows. In Section II,
we describe semi-NMF for GBM-based unmixing.
In Section III, we evaluate the proposed method using

0196-2892 © 2013 IEEE

YOKOYA et al.: NONLINEAR UNMIXING OF HYPERSPECTRAL DATA

1431

synthetic datasets, comparing it with conventional methods.
We present the experimental results for real data in Section IV,
and conclude this paper in Section V.
II. GBM-BASED U NMIXING VIA S EMI -N ONNEGATIVEM ATRIX FACTORIZATION

Halimi et al. [12] proposed three optimization methods
using different algorithms: the Bayesian algorithm using theMarkov-chain Monte Carlo method; the Fan-FCLS algorithm;
and the gradient descent algorithm (GDA). The GDA showed
good results for both synthetic and real hyperspectral datasets.

A. Generalized Bilinear Model

B. Semi-Nonnegative Matrix Factorization for GBM

The BMM takes account of second-order photon interactions between different D endmembers as additional terms in
the LMM, assuming that third or higher order interactions are
negligible [5]. In the BMM, the observed L spectrum of asingle pixel y ∈ R L×1 is given bybi, j ei  e j + n
y = Ea +In this paper, the new optimization method based on semiNMF is introduced to speed up the optimization of a whole
image in matrix form. The observed hyperspectral image canbe reshaped as a matrix form Y ∈ R L×P , with P representing
the number of pixels. The BMM for the whole image is given
in matrix form by
Y = EA + MB + N(i, j )∈G

G = {(i, j )|i < j, i ∈ {1, . . . ,D} , j ∈ {1, . . . ,D}}

(1)

(5)

where A ∈ R D×P is the abundance matrix with each column

where E ∈ R L×D is the endmember matrix with the i thDcolumn vector {ei }i=1
∈ R L×1 representing the i th endmember
D×1
spectrum, a ∈ R
is the abundance vector, bi, j is theinteraction abundance between the i th and j th endmembers, is the Hadamard (element-wise product) operation, andn ∈ R L×1 is the additive noise. On the right-hand side, the firstterm denotes linear mixing, and the second term represents
a second-order interaction, which is assumed to be a linear
combination of the cross-products of the interacting virtualendmembers, i.e., the bilinear spectral mixing. The constraint
i < j defined in (1) is based on the commutative propertyof the Hadamard operation, i.e., ei  e j = e j  ei . On theassumption that an endmember spectrum is a representativeof the multi-order self-interaction, the i th pure endmember
spectrum ei = (ei1 , . . . , ei L )T ∈ R L×1 is defined by a power
series of the self-interactionP
vector {al }l=1
∈ R D×1 representing the abundance vector at
the lth pixel, M ∈ R L×D(D−1)/2 is the bilinear endmember
matrix, B ∈ R D(D−1)/2×P is the interaction abundance matrixP∈ R D(D−1)/2×1 representing
with each column vector {bl }l=1
the interaction abundance vector at the lth pixel, and N ∈R L×P is the noise matrix. GBM-based unmixing becomes thefollowing minimization with respect to A and B:∞


Y1 = EA + N
Y2 = MB + N.

eik =αn fiknminimize Y − EA − MB2F
subject to A  0,ai ≥ 0 i ∈ {1, . . . , D} and0 ≤ ci, j ≤ 1i< ji ∈ {1, . . . , D} ,0BA∗ (7)

i=1(i, j )l = (A)il (A) j l (l ∈ {1, . . . , P}), the operator
· F denotes the Frobenius norm, and the symbol denotes
inequality component-wise. By introducing Y1 = Y−MB andY2 = Y − EA, (5) is written as

n=1

D


(A)il = 1,

(6)where (A∗ )(2)where eik (k ∈ {1, . . . ,L}) is the kth element of theendmember spectrum ei , αn is the abundance of nth-orderself-interaction, and fi = ( f i1 , . . . , f i L )T ∈ R L×1 is thefirst-order self-interaction spectrum of the i th endmember.
Consequently, bi,i corresponding to the second-order selfinteraction is eliminated in (1) because the second-order selfinteraction term, i.e., n = 2, is inherent in the endmember
matrix E.
From a physical perspective, the GBM introduces the nonlinear mixing coefficient bi, j as bi, j = ci, j ai a j and assumes
the following constraints:

D


(8)
(9)

Alternately minimizing Y1 − EA2F and Y2 − MB2F withrespect to A and B, respectively, can be used for the local
optimization of the original problem (6) performed using
an alternating optimization algorithm. Owing to physicalconstraints, all the components of E, M, A, and B are
nonnegative. Therefore, the minimization of (6) can be doneusing alternative update rules of semi-NMF [23]. Semi-NMF
is a problem that factorizes the nonrestricted matrix X into thenonrestricted matrix F and the nonnegative matrix G such that
X = FGT [23]. Ding et al. [23] proposed an iterative updating
algorithm that alternatively updates F and G using−1

ai = 1

(3)

i=1

(10)
F = XG(GT G)

 

+
−
−
+
G ← G· ∗ (XT F) + G(FT F) · / (XT F) + G(FT F)(11)j ∈ {1, . . . , D} .

(4)

When the endmember spectra are known, GBM-based unmixing turns to the optimization of the abundance (a) and theinteraction coefficient (c) under the constraints of (3) and (4).

where · ∗ and · / denote multiplication and division elementwise. (C)+ and (C)− are the positive and negative parts of thematrix C defined as C+ = (|C| + C)/2, C− = (|C| − C)/2.
Semi-NMF optimization is guaranteed to converge to a local
optimum using alternative updating rules.

1432

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 52, NO. 2, FEBRUARY 2014

When E is given and M is calculated from E, GBM-based
unmixing is solved using the following update rules for A
and B:AT ← AT · ∗


 

+
−
−
+
(Y1T E) + AT (ET E) · / (Y1T E) + AT (ET E)

(12)

 

+
−
−
+
(Y2T M) +BT (MT M) · / (Y2T M) +BT (MT M) .

BT ← BT · ∗(13)Since semi-NMF converges to a local optimum, initialization is critical to the final results. We assume that theLMM can be an approximation of the BMM. The fully
constrained least-squares (FCLS) method [24], which is awell-known abundance estimation method based on the LMM,
is used to initialize the abundance matrix A, and the interaction abundance matrix B is initialized as δ × A∗ with asmall δ.
For any observed spectrum (y), the linear mixing term islarger than the bilinear mixing term
Ea − Mb D

a i ei −

a i ei 1 −

i=1

=

D

i=1

a i a j ei  e j(i, j )∈G

i=1
D




a i ei1 + ai2

j =i a j



2
 0.

(14)

When the signal-to-noise ratio (SNR) between Ea + Mb
and n is larger than 2, the SNR between Ea and nis larger than 1. Generally, the SNR of hyperspectral sensors
is much larger than 2; therefore, all components of Y1 can beassumed to be nonnegative. In this case, the minimizing thecost function Y1 − EA2F can be done by NMF [17], which
factorizes a nonnegative matrix into two nonnegative matrices.
A can converge to a local optimum by the NMF update rule
given by A ← A· ∗ ET Y1 · / ET EA .
(15)
If any component of Y1 is negative owing to extreme noise at
the beginning, it should be initialized to zero.
The convergence property based on the update rulesof (12) and (13) or (15) and (13) is proved as follows.
In [23], it is guaranteed that, by fixing E, the residual error
Y1 − EA2F (Y − EA − MB2F ) decreases monotonicallyunder the update rule for A. In the same way, by fixing M, theresidual error Y2 − MB2F (Y − EA − MB2F ) decreases
monotonically under the update rule for B. Therefore,Y − EA − MB2F decreases monotonically under the update
rules (12) and (13). In [18], it is also proved that, when Y1
is nonnegative, Y1 − EA2F (Y − EA − MB2F ) decreases
under the update rule (15). Therefore, Y − EA − MB2F
decreases monotonically with the alternative updating of(15) and (13). The algorithm of GBM-based unmixing usingsemi-NMF is as follows.

Algorithm 1 Generalized Bilinear Model-Based Nonlinear
Unmixing Via Semi-NMF
Input: Hyperspectral data Y ∈ R L×P and endmember matrix
E ∈ R L×D .
Output: Abundance matrix A ∈ R D×P and interaction abundance matrix B ∈ R D(D−1)/2×P .
Step 1: A is initialized using FCLS based on the LMM.
Step 2: A∗ is calculated and B is set as δ × A∗ with a smallδ.
Step 3: A and B are alternatively updated using (12) and (13)(semi-NMF only) or (15) and (13) (NMF and semi-NMF). If
any element of B exceeds that of A∗ , it is replaced with thatof A∗ .
To satisfy the abundance sum-to-1constraint, the methodcited from [24] is adopted.
III. E VALUATION OF S YNTHETIC DATA
Synthetic hyperspectral images (20 ×20 pixel size) are used
to evaluate the unmixing performance of the proposed method,
considering the robustness of the algorithm for the accuracy ofendmember spectra and the collinearity issue. Since bilinear
spectral mixing often occurs between vegetation and soil
[1]–[5], four spectra, i.e., the spectra of two types of plants,
soil, and water in the U. S. Geological Survey spectral library,
are selected as materials of synthetic images. The leaves ofspurge and oak are chosen as vegetation, and the Stonewall
Playa surface is used as soil. The original endmember spectra
are spectrally downsampled with the spectral characteristicof the Airborne Visible/Infrared Imaging Spectrometer. To
evaluate the robustness of the algorithm for spectral complexity and the collinearity issue, two types of synthetic datasetwith different combinations of endmembers are generatedusing three (spurge leaf, soil, and water) and four endmember
spectra, respectively. For each combination of endmembers,
we use three different spectral mixing models, i.e., the LMM(image 1), the GBM (image 2), and a hybrid of the LMM
and GBM (image 3). The half pixels in the hybrid image are
generated by the LMM and GBM. To obtain synthetic images
without pure pixels, all abundance fractions are generated
uniformly on a simplex with a cutoff threshold fixed at 0.8. For
the GBM, the interaction coefficients are set using uniformly
distributed random values in [0 1]. Gaussian noises are added
to all images, setting SNR to 20, which is generally a noise
condition worse than that for real hyperspectral data. Two
scenarios are used for endmember extraction. In the firstscenario, true endmember spectra are given. In the second
scenario, vertex component analysis (VCA) [25], which isone of the most common convex-geometry-based endmember
extraction methods that entail a pure-pixel assumption, is used
for endmember extraction.
First, we compare the performances of the proposed two
methods, i.e., the semi-NMF-only optimization and the NMF
and semi-NMF optimization, using synthetic datasets generated using three endmembers. We adopt two criteria, i.e., thereconstruction error (RE) of unmixing, and the root mean

YOKOYA et al.: NONLINEAR UNMIXING OF HYPERSPECTRAL DATA

1433

NMF & Semi-NMF

Semi-NMF
1st Scenario
x 10

1.124

Image 1

-2

x 10

1.20

-2

Image 2x 10

-2Image 3

1.15

RE

1.1235

1.18

1.123

1.16

1.1225

1.14

1.14

1.13

1.12

1.12

1.122
1.1215
0
4.10

200

x 10

400

0

600

-3

200

x 10

400

600

x 10

200

400

600

200

400

600

-2

1.7

2.5
4.05

0

-2

2.4
1.6

2.3

RMSE
4.00

2.2

1.5

2.1
3.95

2.0
0

200

400

600

1.9

1.4
0

200

400

600

0

Number of itera ons

2nd Scenario
1.42

x 10

-2

Image 1

1.40

RE

1.38
1.36
1.34
0

200
x 10

400

600

-2

Image 2
1.50

1.41

1.45

1.40

1.40

1.39

1.35

1.38

1.30

1.37

0

8.47

200
x 10

x 10

1.55

1.42

-2

8.92

RMSE

x 10

1.43

400

600

1.25

0

-2

7.8

8.46

-2

x 10

Image 3

200

400

600

200

400

600

-2

7.7

8.45

8.90

7.6

8.44
8.88

8.43

7.5

8.42

8.86
0

200

400

600

8.41

7.4
0

200

400

600

0

Number of itera ons
Fig. 1.

Changes in RE of unmixing and RMSE of abundance with number of iterations.

square error (RMSE) of abundance fractions. The RE of
unmixing is calculated to measure the distance between theobserved (yl ) and estimated (yl ) spectra as
1
RE = 
LPPyl − yl 2 .

(16)

l=1

The RMSE of abundance fractions is the most importantcriterion for evaluating unmixing performance, which can beused to compute the accuracy of estimated abundance maps
using the actual (ai ) and estimated (ai ) abundances as
1
RMSE = 
DP

D

i=1

ai − ai 2 .

(17)

Fig. 1 shows the changes in the RE of unmixing and in the
RMSE of abundance fractions with the number of iterations
for three synthetic images. In the first scenario, the NMF
and semi-NMF optimizations show faster convergence forimages 2 and 3 in the RE of unmixing and smaller RMSEs
of abundance fractions, which monotonically decrease. In thecase of image 1, which is based on the LMM, the RMSE
of abundance fractions increases in both methods because
the FCLS method converges to a global optimum. In thesecond scenario, the NMF and semi-NMF optimizations show
faster convergence in the RE of unmixing for all the datasets;
however, the RMSE of abundance fractions increases afterconvergence for image 2. This means that overfitting occurs
in the optimization when the estimated endmember spectra

1434IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 52, NO. 2, FEBRUARY 2014

D=6

D=16

D=11

A

D=21
0.12

B
C

TABLE IC OMPARISON OF U NMIXING P ERFORMANCES W ITH T HREEE NDMEMBERS (V EGETATION , S OIL , AND WATER ) U SING(a)

A CTUAL AND E STIMATED E NDMEMBERS (F IRST AND
S ECOND S CENARIO ): RE S AND SAM S OF U NMIXING

0

AND RMSE S OF A BUNDANCE F RACTIONS

0.2

First Scenario
FCLS

RE (×10−2 )GDA
Semi-NMF

Image 1

1.212

1.212

1.210

0.409

0.410

0.416

Image 2

1.280

1.255

1.216

2.425

2.356

2.169

Image 3

1.253

1.236

1.216

1.691

1.650

1.519

(b)
0

Fig. 2.
(a) RSS maps of second-order interaction terms. (b) Unmixing
residual errors.

RMSE (×10−2 )FCLS
GDA
Semi-NMF

Second Scenario
FCLS

RE (×10−2 )GDA
Semi-NMF

FCLSImage 1

1.417

1.305

1.293

7.480

7.201

7.348

Image 2

1.426

1.378

1.374

8.448

8.388

8.435

Image 3

1.465

1.314

1.298

7.739

7.440

7.408

0.7
tree
tree (shadow)
grass 1
grass 2
grass 3
grass 4
soil 1soil 2road 1road 2

0.6

Reflectance

0.5
0.4
0.3

TABLE II

0.2

C OMPARISON OF U NMIXING P ERFORMANCES W ITH F OUR E NDMEMBERS(T WO V EGETATION , S OIL , AND WATER ) U SING A CTUAL AND E STIMATED

0.1
0
400

500

600

700
800
900
Wavelength (nm)(a)
Fig. 3.

RMSE (×10−2 )GDA
Semi-NMF

1000

1100

(b)

E NDMEMBERS (F IRST AND S ECOND S CENARIO ): RE S AND SAM S OFU NMIXING AND RMSE S OF A BUNDANCE F RACTIONSFirst Scenario(a) Endmember spectra and (b) location of endmembers.

contain errors. Hereafter, we refer to the NMF and semi-NMF
method as semi-NMF. To avoid overfitting, early stoppage isadopted in the optimization by setting the number of iterations
to 300.
Secondly, semi-NMF is compared with the conventional
methods, i.e., FCLS and GDA, which are based on the
LMM and GBM, respectively. To evaluate the robustness forspectral complexity and the collinearity issue, we use two
types of synthetic dataset generated with different numbers ofendmembers. Tables I and II show the results of unmixing
performance for the synthetic datasets with three and four
endmembers, respectively. FCLS shows the smallest RMSE
of abundance fractions for image 1 in the first scenario.
In Tables I and II, semi-NMF shows the smallest REs of
unmixing for all the images in both scenarios and the smallest
RMSEs of abundance fractions for images 2 and 3 in thefirst scenario. This indicates that the semi-NMF method canconverge to a good local optimum with accurate endmember spectra. When endmember spectra contain errors in thesecond scenario, the RMSE of abundance fractions markedly
increases for all the images, and even the REs of unmixing
are comparable with those in the first scenario. For both
combinations of endmembers (Tables I and II), the differencein RMSE between the first and second scenarios is much larger
than the difference in RMSE between the linear and bilinear
spectral unmixing methods. This suggests that the accuracy ofendmember spectra has a larger effect on the final abundance
estimation than the second-order scattering effect. When thenumber of endmembers increases, as shown in Table II, in thefirst scenario, the REMS of abundance fractions increases for

FCLS

RE (×10−2 )GDA
Semi-NMF

RMSE (×10−2 )FCLS
GDA
Semi-NMF

Image 1

1.449

1.449

1.447

1.985

1.985

2.002

Image 2

1.483

1.494

1.450

6.053

5.902

5.390

Image 3

1.466

1.471

1.452

4.449

4.342

4.035

Second Scenario
FCLS

RE (×10−2 )GDA
Semi-NMFFCLSRMSE (×10−2 )GDA
Semi-NMF

Image 1

1.544

1.569

1.480

7.325

7.285

7.258

Image 2

1.499

1.520

1.427

7.627

7.487

7.484

Image 3

1.549

1.571

1.429

8.385

8.187

8.151

all the methods, which may be due to the collinearity issuecaused by the highly correlated spectra between endmembers
and interacting virtual endmembers. Even in this case, withaccurate endmember spectra, GDA and semi-NMF outperform
FCLS for images 2 and 3. This confirms that GBM-based
unmixing has potential for dealing with spectral complexity.
According to the smallest RMSEs of abundance fractions,semi-NMF is a promising optimization method for GBMbased nonlinear unmixing.
IV. E XPERIMENT FOR R EAL DATA
GBM-based unmixing using semi-NMF is applied to airborne hyperspectral data to evaluate the numerical unmixing performance for real data, to examine the impact ofendmember extraction on bilinear spectral unmixing, and todemonstrate the difference in the abundance map obtained
between the LMM and GBM. The dataset was acquired
using the compact airborne spectrographic imager (CASI-3)
on June 19, 2009. The study area was the agricultural field in
the Motonopporo farm in Hokkaido, Japan, because multiple

YOKOYA et al.: NONLINEAR UNMIXING OF HYPERSPECTRAL DATA

1435

TABLE IIIC OMPARISON OF RE S (×10−3 ) AND SAM S ( IN D EGREE ), AND C OMPUTATIONAL C OSTS (IN S )

(a)

RE

D=6
SAM

Time

RE

D = 11
SAM

Time

RE

D = 16
SAM

Time

RE

D = 21
SAM

Time

FCLS

10.56

2.695

20.5

6.320

1.824

57.4

5.025

1.488

95.2

4.738

1.433

143.5

GDA

10.55

2.692

122.2

6.318

1.823

223.1

5.048

1.479

670.6

4.721

1.431

785.7

Semi-NMF

10.26

2.640

55.4

6.258

1.816

106.6

4.975

1.483

169.7

4.701

1.429

250.4

tree

tree (shadow)

grass 1

grass 2

grass 3

1

grass 4

soil 1soil 2road 1road 2
0

(b)

tree

tree (shadow)

grass 1

grass 2

grass 3-2

×10

2

grass 4

soil 1soil 2road 1road 2-2

Fig. 4.

(a) Abundance maps and (b) differences of abundance maps with LMM.

scattering is a common phenomenon between vegetation andsoil [1]–[5]. The area consisted of a devastated grassland,
trees, unpaved rural roads, an agricultural field, and an asphalt
road. We selected a 150×150-pixel-size image with 68 spectral
channels over 410–1070 nm and a 1-m ground sampling
distance. The data initially measured as radiance was converted
into reflectance.
The performance of the proposed method is numericallyevaluated by the RE and the spectral angle mapper (SAM)of unmixing, and the computational time of Matlab implementations was compared with that of FCLS and GDA. The
CPU used is Intel Core i7 CPU 2.80 GHz, with a memory
capacity of 16 GB. SAM is calculated using the observed (yl )and estimated (yl ) spectra as
1 yl · ylarccos.
yPl  yl 
P

SAM =

l=1

(18)

RE and SAM are common criteria for evaluating theperformance of the unmixing procedure for real datasets.
By changing the number of endmembers extracted by VCA
and comparing the second-order interaction term and theresidual errors of the GBM, we examined the impact ofendmember extraction on bilinear unmixing. The effects of
the interaction term and the residual errors of unmixing ata single pixel are calculated by the root-sum-square (RSS)
method. After the endmember extraction, abundance maps andthe impact of bilinear spectral mixing on these abundancemaps are calculated.
Table III shows a comparison of the RE and SAM values,
and the computational cost determined using four different numbers of endmembers, i.e., D = 6, 11, 16, and 21.
The Semi-NMF method shows the best performance in termsof RE and SAM for many different numbers of endmembers.
In addition, the computational time is smaller than for GDA.
In particular, even when the number of endmembers increases,

1436

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 52, NO. 2, FEBRUARY 2014

the semi-NMF method converges to good local optima at
reasonable computational cost owing to the simple update
rules in matrix form. This result indicates that semi-NMF isa practical method of solving the optimization of GBM-based
unmixing owing to its convergence performance and simple
implementation.
Fig. 2 shows RSS maps of the second-order interaction termand unmixing residual errors. As the number of endmembers
increases (6, 11, 16, and 21), the two maps converge. We concluded that 16 endmembers are sufficient. In spite of the large
number of endmembers, the bilinear spectral mixing effectclearly appears in areas B and C in Fig. 2, i.e., the boundary
between grass and trees, and the narrow path, where there
is a 3-D multilayered structure, which is consistent with thephysical model of GBM-based unmixing. With D= 6, area
A in Fig. 2, corresponding to grass, shows a large bilinear
spectral mixing effect; however, the residual error still remains
large in this area. It significantly decreases after grass is
detected as the endmember at D = 11. This means that
endmember extraction has a larger influence on residual errorsof unmixing than the bilinear mixing effect, which is also
confirmed from the difference between the first and second
scenarios in the simulation of synthetic datasets. Therefore,
to discuss nonlinear spectral mixing, accurate endmember
extraction is necessary.
Since VCA is sensitive to the outlier, it can detect someendmembers that are only present in very specific regions.
To visualize the abundance map and the bilinear mixing effect,
we manually eliminated such type of specific endmembers andreprocessed GBM-based unmixing with 10 endmembers, as
shown in Fig. 3(a). In Fig. 3(b), the colored circles indicate thelocations of these endmembers, with each color corresponding to that of the spectrum in Fig. 3(a). The endmembers
are labeled by visual judgment considering the locations ofendmembers and comparing the abundance maps and RGB
color image. Fig. 4(a) and (b) show the abundance maps forthese 10 materials and their differences from the LMM, which
demonstrate the impacts of the second-order interaction onabundance maps. In some regions, the abundance fractions
changed by more than 2% compared to the LMM becauseof the multiple scattering effects. Grass 4 and soil 2 show
relatively large changes from the LMM at their boundary
owing to the second-order scattering interaction, which is consistent with previous research reports on multiple interactions
that often occur between vegetation and soil. This qualitative
result implies that GBM-based unmixing using semi-NMF has
potential for dealing with the spectrally complex BMM in
hyperspectral unmixing.
V. C ONCLUSION
In this paper, we presented a novel optimization method
based on semi-NMF for GBM-based nonlinear unmixing.
Semi-NMF enables the optimization of GBM-based unmixing
to process a whole image in matrix form with simple update
rules. The performance of the proposed method is evaluatedusing synthetic hyperspectral datasets, considering the robustness for the error of endmember extraction and the collinearity

issue, and compared with the LMM-based method and theconventional optimization for GBM-based unmixing. Withaccurate endmember spectra, GBM-based unmixing showedrobust results for the spectral complexity and collinearity issue.
In particular, the semi-NMF method gave the smallest RMSE
of abundance fractions. The errors of endmember extraction
resulted in RMSEs of abundance fractions being larger than the
difference in RMSE between the LMM and the GBM. Therefore, accurate endmember extraction is necessary to discuss
nonlinear unmixing. We applied GBM-based unmixing usingsemi-NMF to a real airborne hyperspectral image taken overan agricultural field with many endmembers. The semi-NMF
method could be used to visualize the interaction between
vegetation and soil and its effect on final abundance maps
at reasonable computational cost.
ACKNOWLEDGMENT
The authors would like to thank the reviewers for the many
valuable comments and suggestions. N. Yokoya acknowledges
the support from the Japan Society for the Promotion of
Science.
R EFERENCES[1] D. A. Roberts, M. O. Smith, and J. B. Adams, “Green vegetation,
nonphotosynthetic vegetation, and soils in AVIRIS data, “ Remote Sens.
Environ., vol. 44, nos. 2–3, pp. 117–126, May 1993.
[2] C. C. Borel and S. A. W. Gerstl, “Nonlinear spectral mixing models
for vegetative and soil surfaces,” Remote Sens. Environ., vol. 47, no. 3,
pp. 403–416, Mar. 1994.
[3] T. W. Ray and B. C. Murray, “Nonlinear spectral mixing in desert
vegetation,” Remote Sens. Environ., vol. 55, no. 1, pp. 59–64, Jan. 1996.
[4] L. Zhang, D. Li, Q. Tong, and L. Zheng, “Study of the spectral mixture
model of soil and vegetation in PoYang Lake area, China,” Int. J. Remote
Sens., vol. 19, no. 11, pp. 2077–2084, 1998.
[5] B. Somers, K. Cools, S. Delalieux, J. Stuckens, D. V. Der Zande,
W. W. Verstraeten, and P. Coppin, “Nonlinear hyperspectral mixture
analysis for tree cover estimates in orchards,” Remote Sens. Environ.,
vol. 113, no. 6, pp. 1183–1193, Jun. 2009.
[6] J. F. Mustard, L. Li, and G. He, “Compositional gradients across marehighland contacts: The importance and geological implications of lateralmixing,” J. Geophys. Res.—Planets, vol. 103, pp. 419–425, Aug. 1998.
[7] R. Close, P. Gader, A. Zare, J. Wilson, and D. Dranishnikov, “Endmember extraction using the physics-based multi-mixture pixel model,” Proc.
SPIE, vol. 8515, p. 85150L, Aug. 2012.
[8] R. Close, “Endmember and proportion estimation using physics-based
macroscopic and microscopic mixture models,” Ph.D. dissertation, Dept.
Comput. Eng., Florida Univ., Gainesville, FL, USA, Dec. 2011.
[9] W. Fan, B. Hu, J. Miller, and M. Li, “Comparative study
between a new nonlinear model and common linear model for analysing
laboratory simulated-forest hyperspectral data,” Int. J. Remote Sens.,
vol. 30, no. 11, pp. 2951–2962, Jun. 2009.
[10] J. M. Nascimento and J. M. Bioucas-Dias, “Nonlinear mixture model for
hyperspectral unmixing,” Proc. SPIE, vol. 7477, p. 74770I, Sep. 2009.
[11] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, “Nonlinear
unmixing of hyperspectral images using a generalized bilinear model,”IEEE Trans. Geosci. Remote Sens., vol. 49, no. 11, pp. 4153–4162,
Nov. 2011.
[12] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, “Unmixing
hyperspectral images using the generalized bilinear model,” in Proc.
IEEE Int. Geosci. Remote Sens. Symp., Jul. 2011, pp. 1886–1889.
[13] J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader,and J. Chanussot, “Hyperspectral unmixing overview: Geometrical,
statistical, and sparse regression-based approaches,” IEEE J. Sel. Topics
Appl. Earth Observ. Remote Sens., vol. 5, no. 2, pp. 354–379, Apr. 2012.
[14] B. Hapke, Theory of Reflectance and Emittance Spectroscopy. Cambridge, U.K.: Cambridge Univ. Press, 1993.
[15] S. Weisberg, Applied Linear Regression. New York, NY, USA: Wiley,
1985.

YOKOYA et al.: NONLINEAR UNMIXING OF HYPERSPECTRAL DATA

[16] X. Chen, J. Chen, X. Jia, B. Somers, J. Wu, and P. Coppin, “A quantitative analysis of virtual endmembers’ increased impact on the collinearity
effect in spectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 49,
no. 8, pp. 2945–2956, Aug. 2011.
[17] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative
matrix factorization,” Nature, vol. 401, pp. 788–791, Oct. 1999.
[18] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix
factorization,” in Proc. Conf. Adv. Neural Inf. Process. Syst., vol. 13,
2001, pp. 556–562.
[19] L. Miao and H. Qi, “Endmember extraction from highly mixed datausing minimum volume constrained nonnegative matrix factorization,”
IEEE Trans. Geosci. Remote Sens., vol. 45, no. 3, pp. 765–777,
Mar. 2007.
[20] S. Jia and Y. Qian, “Constrained nonnegative matrix factorization for
hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 47,
no. 1, pp. 161–173, Jan. 2009.
[21] A. Huck, M. Guillaume, and J. Blanc-Talon, “Minimum dispersion
constrained nonnegative matrix factorization to unmix hyperspectral
data,” IEEE Trans. Geosci. Remote Sens., vol. 48, no. 6, pp. 2590–2602,
Jun. 2010.
[22] X. Liu, W. Xia, B. Wang, and L. Zhang, “An approach based on
constrained nonnegative matrix factorization to unmix hyperspectral
data,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 2, pp. 757–772,
Feb. 2011.
[23] C. Ding, T. Li, and M. I. Jordan, “Convex and semi-nonnegative matrix
factorization,” IEEE Trans. Patten Anal. Mach. Intell., vol. 32, no. 1,
pp. 45–55, Jan. 2010.
[24] 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, Mar. 2001.
[25] J. M. Nascimento and J. M. Bioucas-Dias, “Vertex component analysis:
A fast algorithm to unmix hyperspectral data,” IEEE Trans. Geosci.
Remote Sens., vol. 43, no. 4, pp. 898–910, Apr. 2005.

Naoto Yokoya (S’10) received the M.Sc. and Ph.D.
degrees in aerospace engineering from the University of Tokyo, Tokyo, Japan, in 2010 and 2013,
respectively.
He has been a Research Fellow with the Japan
Society for the Promotion of Science, since 2012.
His current research interests include hyperspectral
data analysis.

1437

Jocelyn Chanussot (M’04–SM’04–F’12) received
the M.Sc. degree in electrical engineering from theGrenoble Institute of Technology (Grenoble INP),
Grenoble, France, and the Ph.D. degree from Savoie
University, Annecy, France, in 1995 and 1998,
respectively.
He was with the Geography Imagery Perception Laboratory, Delegation Generale de l’Armement(DGA - French National Defense Department), Bagneux, France, in 1999. Since 1999, he has been with
Grenoble INP, where he was an Assistant Professor
from 1999 to 2005, an Associate Professor from 2005 to 2007, and is currently
a Professor of signal and image processing. He is currently conducting his
research at the Grenoble Images Speech Signals and Automatics Laboratory(GIPSA-Lab). His current research interests include image analysis, multicomponent image processing, nonlinear filtering, and data fusion in remotesensing.
Dr. Chanussot is the founding President of IEEE Geoscience and Remote
Sensing French chapter (2007–2010) which received the 2010 IEEE GRSS Chapter Excellence Award “for excellence as a Geoscience and Remote
Sensing Society chapter demonstrated by exemplary activities during 2009.”
He was the recipient of the NORSIG 2006 Best Student Paper Award, theIEEE GRSS 2011 Symposium Best Paper Award and of the IEEE GRSS 2012
Transactions Prize Paper Award. He was a member of the IEEE Geoscienceand Remote Sensing Society AdCom (2009–2010), in charge of membershipdevelopment. He was the General Chair of the first IEEE GRSS Workshop
on Hyperspectral Image and Signal Processing, Evolution in Remote sensing(WHISPERS). He is the Chair (2009–2011) and was the Cochair of the
GRS Data Fusion Technical Committee (2005–2008). He was a memberof the Machine Learning for Signal Processing Technical Committee of theIEEE Signal Processing Society (2006–2008) and the Program Chair of theIEEE International Workshop on Machine Learning for Signal Processing,(2009). He was an Associate Editor for the IEEE G EOSCIENCE AND R EMOTES ENSING L ETTERS (2005–2007) and for Pattern Recognition (2006–2008).
Since 2007, he has been an Associate Editor for the IEEE T RANSACTIONS ON
G EOSCIENCE AND R EMOTE S ENSING. Since 2011, he has been the Editorin-Chief of the IEEE J OURNAL OF S ELECTED T OPICS IN A PPLIED E ARTH
O BSERVATIONS AND R EMOTE S ENSING.

Akira Iwasaki received the M.Sc. degree inaerospace engineering and the Doctoral degree inengineering from the University of Tokyo, Tokyo,
Japan, in 1987 and 1996, respectively.
He joined the Electrotechnical Laboratory, where
he engaged in research on space technology andremote sensing system in 1987. He is currently aProfessor with the University of Tokyo.