2746 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 50, NO. 7, JULY 2012
for an abundance estimation process such as the one conducted
under the assumption that pure pixels are present in the input
data. Recently, a similar approach has been developed using se-
quential linear programming to solve the MVES problem [40],
[44]. As will be shown, standard available nonlinear optimi-
zation algorithms can be used to solve the same problem.
The remainder of this paper is organized as follows.
Section II describes the method and its underlying principles.
Section III provides a detailed experimental assessment of the
performance of the method using a database of synthetic hyper-
spectral scenes generated using fractals and a real hyperspectral
data set collected by AVIRIS. In all cases, the performance is
compared to that of state-of-the-art techniques for endmember
extraction which assume the presence of pure pixels [OSP,
N-FINDR, VCA, spatial–spectral endmember extraction
(SSEE), AMEE, SPP, and RBSPP] and also to that of
techniques that do not incorporate such assumption (MVES
and MVSA). Section IV concludes with some remarks and
hints at plausible future research.
II. Minimum-Volume Simplicial Enclosure Method
A. Principle of Volume Reduction
Linear spectral unmixing methods [3] exploit the idea that
the observations x in m-dimensional space (m = 224 spectral
bands for the AVIRIS imaging spectrometer) are in reality com-
binations of only n endmembers. This means that observations
in the space E = {x = Ea|a ∈ R
n
,
i
a
i
=1}⊂R
m
can be
represented by points in V = {z = Va|a ∈ R
n
,
i
a
i
=1}⊂
R
n−1
, which is much lower in dimension, i.e., n m, where
V =[v
1
,...,v
n
] is an (n − 1) × n matrix of the so-called
vertices v
i
. The transformation is based on some observations
from linear algebra and works as follows.
Let E
−
be the matrix E
−
=[e
2
− e
1
,e
3
− e
1
,...,e
n
− e
1
].
Its columns span the linear space E
−
such that one can
write E as E = {x = e
1
+ E
−
y|y ∈ R
n−1
}. More precisely, let
Γ be a matrix where its columns form an orthonormal basis
of E
−
. The number of columns of Γ is equal to the rank
of E
−
. The rank is n − 1 if the columns of E are affine
independent. This means that the space E can be described as
E := {x = b +Γd|d ∈ R
n−1
}, with b as an arbitrary given m-
vector in E and Γ as an arbitrary m × (n − 1) orthonormal
matrix, such that Γ = E
−
, i.e., its columns have length 1
and are orthogonal with respect to each other. Now, given a
choice for b and Γ, let the vertices v
i
be solutions of
Γv
i
= e
i
− b, i =1,...,n. (2)
Then, in matrix notation, we have E =ΓV + b1
T
, where 1
is the all-one vector of appropriate dimension. For each abun-
dance vector a (keep in mind that its elements sum to one),
we have a unique point x ∈ E as well as z ∈ V related as
x = Ea =(ΓV + b1
T
)a =Γz + b.
Two questions arise at this point: 1) What choices are conve-
nient for b and Γ? and 2) what are the consequences for model
(1), which includes white noise, if we have such a vector b
and orthonormal basis matrix Γ? To start with the latter, let us
consider the so-called score vector z in (n − 1)-dimensional
space
z = Va+ ξ (3)
where Γξ is the projection of on the space Γ. It can be de-
rived that ξ is white noise. Due to Γ being orthonormal, Γ
T
Γ is
the unit matrix, and the projection leads to ξ =(Γ
T
Γ)
−1
Γ
T
=
Γ
T
. Notice that the variance–covariance matrix of ξ is Γ
T
Γ,
i.e., the unit matrix. The relation of (3) with the original model
(1) is
(x − b)=Ea − b + =ΓVa+
=Γ(Va+ ξ)+ζ =Γz + ζ (4)
where ζ is the projection of on the orthoplement of Γ.So
in fact, is decomposed into =Γξ + ζ. The consequence
of these algebraic observations and the underlying stochastic
model is that, given a support vector b and matrix Γ, the original
model (1) in m-dimensional space reduces to model (3) in
(n − 1)-dimensional space. The observations X =[x
1
,...,x
r
]
that, without noise, would lay in the simplex constituted by the
so-called convex hull of E, conv(E) ⊂ E, have a one-to-one
relation with their representation Z =[z
1
,...,z
r
] in simplex
S =conv(V ). Given these observations, a desirable objective
is to find out the endmember representation V and the original
endmember signals in E; this process is called endmember
identification.
In [41], a matrix factorization approach is used that simulta-
neously minimizes least squares and the simplex volume for
endmember generation. For volume reduction, it is good to
add a direction without variation to t he principal components;
zero variation gives zero volume. In turn, the least squares
idea looks for directions with maximum variation. Miao and
Qi [41] use weights to deal with these two objectives. In
this paper, we follow a similar strategy based on estimating
the number of interior/exterior points by weighting these two
terms. Specifically, we investigate a procedure which takes a
hierarchical approach by first estimating the subspace E in
which the n endmembers are lying by means of PCA ( a similar
approach was adopted in [45]).
The choice of using PCA means that space E = b + Γ
is estimated by taking b =
x =(1/r)
j
x
j
and Γ=C =
[c
1
,...,c
n−1
] is an m × (n − 1) matrix of principal compo-
nents; we have the biggest variation in direction c
1
, the second
biggest in direction c
2
, etc. Following this procedure, the (un-
known) scores z
k
are estimated from observations x
k
by taking
z
k
=(C
T
C)
−1
C
T
(x
k
− x)=C
T
(x
k
− x). Then, the method
that we investigate minimizes (in that subspace) the volume of
the resulting simplex such that it encloses the scores z
k
of the
observed bands of the pixels x
k
.
B. Principle of MVES
The problem of finding the MVES of a set of points z
k
,k =
1,...,r in (n − 1)-dimensional space can be written as a so-
called nonlinear optimization problem [46] with (n − 1) × n