Remote Sensing Analyis of Alteration Minerology Associated With Natural Acid
in the Grizzly Peak Caldera, Sawatch Range, Colorado
David W. Coulter
Ph.D. Thesis Colorado School of Mines, 2006
West Red (Ruby Mountain) from Enterprise Peak
West Red Fe Endmembers from Aviris
Red: Hematite Green: Goethite Blue: Jarosite
(low pH) (high pH)
2 Cross-Correlation Spectral Matching
Cross-Correlation Spectral Matching
wavelength vector with
unknown (observed) image spectrum at given pixel
known reference spectrum; (to be matched against y)
Linear Correlation Coefficient (Pearson’s ):
. Generalize: cross-correlation at position
Measures of fit:
relative to 1 is not particularly good, but
is distributed in null case as Student’s
For non-uniform weights, weight the sums in eq. (1) by
where the incomplete gamma function
If e.g. ,
then the goodness-of-fit is believable. If
is larger than, say, 0.001, the fit may be acceptable if the the
errors are non-normal or have been moderately underestimated. If
then question the model and/or estimation proceedure eq. (2). If latter, consider
Other useful measures Skewness
¡1-¿ Skewness: cross-correlogram of perfect match is parabola:
where is auto-correlogram
of z with itself, is number
of match positions,
is match number.
Another possibly useful statistical measure, if
is binormal and
large , is Fisher’s
Then if you measure the correlation coefficients of the unknown signal y against two reference
the significance of a difference between the two measured correlation coefficients
is given by the complementary error function [6, eq. (13.7.10)]
This can be useful when trying to assess whether a parameterized fit changes
”significantly” when a given change of parameter produces the two references
Non-parametric or Rank Correlation (Spearman’s Rank-Order Correlation
Coefficient). See [6, Sec. 13.8] and .
3 Optimized Cross-Correlation Spectral Matching
Optimized Cross-Correlation Spectral Matching In ”traditional” CCSM, the reference spectrum
z to which
is compared is taken to be a single (pure) endmember. Next assume M endmembers
of interest, and linear mixing. Seek weighting factors
(synthesized pixel intensity at band )
that maximize the cross correlogram
at match position , places
that maximum at , and
minimizes its skew, subject to
and . Significance of
different values of
resulting from different choices of endmember set
be assessed using eq. (10). Coulter  chooses endmember set that maximizes
written (12) is the cross-correlation of two unit vectors:
which is independent of the normalization of
is spectrum at a particular image pixel, then the unit normalization
somewhat independent of shadow and striping.
Equation (14) makes
non-linear in the ; eq (13)
may be maximized w.r.t the
by the non-linear constrained optimizer of one’s choice. It does not
give much insight into the relative band-to-band errors inherent in
written, it assumes they are all equal, and we can pretend:
Given spectrometer calibration, we can do better than this.
4 Linear and Non-Linear Least-Squares Fitting
Least-Squares Fitting Since
is expressed in terms of unit vectors, maximizing
equivalent to minimizing
Minimizing the squared-difference of unit vectors
in (16) is still a
non-linear problem in the
because of how they appear in the denominator of the normalization of
(14). However, if we relax the restrictions that z be a unit vector (and that
can define the equivalent problem
which is linear in the ,
but not yet quite as general as we want. At each band
(at least) three uncertainties in the sensor measurement: those in the received radiance
, those in the center
and those in the sensor’s bandwidth, assumed as a fwhm of a gaussian. Fwhm is
typically on order of band-to-band wavelength spacing. Fwhm errors are usually of
secondary significance and will be ignored.
If we allow for uncertain band centers
and assume normal independent distributions, we can write
as joint probability sensor records band radiance
at recorded wavelength
, when it actually received
(unknown) radiance at actual
(and unknown) wavelength .
Joint probability of pixel’s measured spectrum across N bands is
Associate the (unknown) actual spectrum
with the modeled
represents a constrained optimization problem wherein we wish to minimize
as a function
. The motivation for including
the is to (hopefully) obtain
in the optimal solution.
Eqs (28) and (30) make eq. (27) quadratic in the
the optimization including center wavelengths uncertainties is non-linear. The
can be computed
efficiently enough, but the
will need to be reconvolved against the new center wavelengths for each new set of center
required during the optimization. The problem will obviously be much
simpler if we assume the provided center wavelength values are
”good enough” and simply ignore eqs (27). Dropping the superscript
equations (26) then become
which are the usual linear least-squares normal equations for the
what its worth, for the unconstrained problem the minimal solution vector
will be unique
provided the matrix
5 Image Classification and Statistics
5.1 Image Mean Vector and Covariance Matrix
Image Mean Vector and Covariance Matrix [11, White 2005] If image has
pixels, the mean vector is
is the jth pixel vector of the image
The image covariance matrix is
where is the
matrix of pixel
vectors each of length
is the matrix
mean vectors (
matrix of ones.
5.2 Principal Component Transformation
Principal Component Transformation [7, Smith 1985] Karhunen-Loeve Transformation[12,
White 2005] GRASS imagery i.pca: Let
input-image multi-pixel vector (
mean vector matrix,
output-image multi-pixel vector,
matrix whose rows are the eigenvectors of the covariance matrix ,
arranged by decreasing magnitude of eigenvalue, as typically returned by
the symmetric positive-definite image covariance matrix
are its orthonormal eigenvectors with eigenvalue .
impose ordering on transformed component vectors
. Those with
the largest s.t.
are the PrincipalComponents.
should be related to the noise floor.
5.3 Minimum Noise Fraction
Minimum Noise Fraction[10, pg. 38][3, ] We wish to find a particular coefficient matrix
in some sense maximizes the image S/N, assuming the image pixel vectors
the sum of uncorrelated signal and noise:
where is generalized
eigenvalue of wrt
corresponding generalized eigenvectors. Compare with PCA:
Green suggests be
of unit variance and band-to-band uncorrelated.
Divide image into small
pixel blocks (4x4, 5x5,...)
For each (block, band) get local mean and variance:
is the local
into classes between band min and max values.
bin with most blocks represents mean noise of image. Hope this bin is same for all
Local Means and Local Variances (con’t)
Suppose ”most popular” bin is a P-cube, each side ranging
points. Then the average value over the bin
is desired noise covariance matrix.
Caveat: Assumes image is ”slowly varying enough” that enough of the
blocks are homogeneous, i.e. their covariance is really due to noise, not true
features in the image. Blocks and bins must both be ”small enough” – but not too
Other methods: ”unsupervised training” derived-endmember classification schemes e.g.
LAS’ search ([13, ]) and GRASS’ cluster/maxlik are based upon local
6 Vertex Component Analysis
Vertex Component Analysis is an “unsupervised training” derived-endmember
classification scheme. We use the notation of [5, Nascimento and Dias], which is
slightly different than that used in previous sections.
6.1 VCA Algorithm
Assume linear mixing, and let
number of pixels in image
number of spectral bands recorded at each pixel
number of endmembers, ,
for hyperspectral images.
is an endmember spectral vector of length ,
is a scale factor modeling illumination variability due to surface topography
is the abundance vector containing the fractions of each endmember at pixel
is the “true” noise-free spectral vector of length
at pixel .
is an Lxp mixing matrix that maps an abundance vector
to a “true” spectral vector .
Then the recorded spectral vector at pixel
may be given by
Our goal is to find a abundance vectors
corresponding to some
endmember set . An
appropriate endmember set
is to be determined as part of the VCA algorithm. Endmember identification – the
matching one or more of the VCA generated endmembers to actual specimen spectral
samples such as from USGS spectral library or field ground-truth sampling – may be
done in a subsequent processing step.
Since the set is a simplex,
the set is also a simplex.
However, even assuming ,
th observed vector set belongs to convex cone
owing to different
scale factors .
But the projective projection of the convec cone
properly chosen hyperplane is a simplex with vertices corresponding to the vertices of simplex
. The simplex
is the projective projection
of the convex cone onto
the hyperplane , where the
choice of assures there
are no observed vectors
orthogonal to the hyperplane: .
VCA algorithm accuracy is dependent on image SNR. Several ways to estimate:
Classic definition: SNR = Power(signal)/Power(noise). Select homogeneous
blocks from image, areas where the imageis expected to be flat so that
variations are due to noise. Calculation of data SNR using a (Mean)/(Standard
Deviation) method for a homogeneous target.
is the average intensity of the ith spectral band in the blobk with
is the standard deviation. Then .
We may manually select homogeneous regions from the image, such as still
water (lakes, ponds) or snowfields, or we may use GRASS i.cluster/i.maklik
to (perhaps) find such regions automatically.
. VCA does an image SVD or PCA as an initial dimensionality reduction
step. We might try looking at the ratio of the Power in the SVD (and/or
PCA) components with large singular values – the ones of interest – with
the Power in the components with small singular values. Compare this with
the ratio of the sum of the large eigenvalues with the sum of the small
eignevalues – or sums of squares of each, to see if there’s an easy way to
compute an “effective SNR” that is usable to select between SVD and PCA
in the first step of the VCA algorithm.
Compare the above with sensor data e.g. AVRIS .rcc files Aviris has very
high SNR, over 500 on a good day[4, Sec. 3.3] Satellites e.g. Hyperion
might have a tenth this, SNR is an important consideration and worth some
synthetic image modelling. Nascimento and Dias use
and SNR estimator. We note with an SNR greater than 500 on a good day,
AVIRIS is considered extemely good.[4, Kruse JPL 2002] But with 224
bands it will spectacularly fail this particular test, as .
We speculate this test is insufficent for any hyperspectral image, and something
better related to the number of significant principal components (usually
d) might be in order. Nascimento and Dias have quite a bit on the topics
of estimating SNR and the virtual dimensionality (VD) of an image. See
[5, eqs 10 - 13] and C.-I Chang and Q. Du, “Estimation of number of
spectrally distinct signal sources in hyperspectral imagery”, IEEE Trans.Geosci. Remote Sens., vol. 42, no. 3, pp. 608-619, Mar. 2004. [1, -I Chang
6.1.2 Algorithm 1: Vertex Component Analyis (VCA)
is the image matrix of
each of length
desired dimensionality of problem, the number of endmembers.
may be estimated as the number of “large” eigenvalues of the image covariance
is the fixed threshold value
of the SNR measured wrt the signal subspace. Different values might be
used for different SNR estimators. See above.
estimated mixing matrix returned by VCA.
is the jth column of
are the ith through kth columns of
is the matrix formed from the columns of
listed in index vector
denotes the pseudo-inverse
of pxp matrix as
computed by SVD.
is what you think it is. Write the SVD as
is the pxp diagonal matrix whose diagonal elements are the (ordered) singular values:
Then the pseudo-inverse
may be written where
and the pseudo-inverse
of any scalar
is the orthogonal projector onto the range of
is the orthogonal projector onto the range of
is the orthogonal projector onto the kernel (null-space) of
is the orthogonal projector onto the kernel (null-space) of
has d non-zero diagonal entries, .
is a pxp block-diagonal matrix with a non-zero dxd upper-left block,
and is zero everywhere else.
See http://en.wikipedia.org/wiki/Moore-Penrose_pseudoinverse In particular “ the
pseudoinverse for matrices related to can be computed by applying the
ShermanMorrisonWoodbury formula to update the inverse of the correlation
matrix, which may need less work. In particular, if the related matrix differs from
the original one by only a changed, added or deleted row or column, incremental
algorithms exist that exploit the relationship.” – a fact that might be quite
is an Lxd matrix whose columns are the first d eigenvectors of either
(which is LxL) or the LxL spectral covariance matrix
is the left eigenmatrix computed by SVD whose columns are ordered (left to right) in
order of decreasing singular value. Its first d columns are the most significant
(principal) relative to the rest. We usually select d such that the singular values
are in some sense “too small to be significant”.
1: Compute SNR 2: if
end if 14: initialize
, a pxp auxiliary
by SVD 15: for ()
22: end for 23: if
then 24: is
estimated mixing matrix. 25: else 26: is
estimated mixing matrix. 27: end if
C.-I Chang and Q. Du. Estimation of number of spectrally distinct signal
sources in hyperspectral imagery. IEEE Transactions on Geoscience andRemote Sensing, 42(3):608–619, March 2004.
David W. Coulter. Remote Sensing Analysis of Alteration MineralogyAssociated with Natural Acid Drainage in the Grizzly Peak Caldera,Sawatch Range, Colorado. PhD thesis, Colorado School of Mines, Golden,
A.A. Green, M. Berman, P. Switzer, and M.D. Graig. A transformation
for ordering multispectral data in terms of image quality with implications
for noise removal. Journal of Geophysical Reseach, 90:797 – 804, 1988.
Fred A. Kruse. Comparison of aviris and hyperion for hyperspectral
José M. P. Nascimento and José M. Bioucas Dias. Vertex component
analysis: A fast algorithm to unmix hyperspectral data. IEEE Transactionson Geoscience and Remote Sensing, 43(4), April 2005.
William H. Press, Brian P. Flannery, Saul A. Teukolsky,
and William T. Vetterling. Numerical Recipes in C. Cambridge University
Press, Cambridge, New York, Port Chester, Melbourne, Sidney, 1988.
M.O. Smith, P.E. Johnson, and J.B. Adams. Quantitative determination
of mineral types and abundances from reflectance spectra using principal
component analysis. IEEE Transactions on Geoscience and Remote Sensing,
36:65 – 74, 1985.
Frank D. van der Meer. Extraction of mineral absorbtion features
from high-spectral resolution data using non-parameteric geostatistical
techniques. International Journal of Remote Sensing, 15:2193–2214, 1994.
Frank D. van der Meer and Steven M. de Jong. Imaging Spectroscopy.
Kluwer Academic Publishers, Dordrecht, Boston, London, 2001.
Frank D. van der Meer, Steven M. de Jong, and W. Bakker. ImagingSpectroscopy: Basic analytical techniques, pages 17–62. Kluwer Academic
Publishers, Dordrecht, Boston, London, 2001.
R. A. White. Image mean and covariance:
R. A. White. Karhunen-loeve transformation:
R. A. White. Search unsupervised training site selection: