Contents

1 Introduction

The Gizzly Peak Caldera is located Colorado’s Mineral Belt approximately 15 miles southeast of Aspen. It is the product of a volcanic eruption $\sim 35$ million years ago that ejected some 600 ${\text{ km}}^3$ of magma. With a volcanic explosive index of 7, it was at least four times larger than Tambora.[2, sec 1.5]

Figure 1: West Red (Ruby Mountain) From Enterprise Peak

Coulter sought to identify acidic conditions by the different weathering states of iron oxide, which forms Hematite (red) in acidic conditions, and Goethite (lite yellow) and Jarosite (dark yellow) as pH increases, Jarosite being the most immediate oxidation product of iron sulfide (pyrite).

Figure 2: West Red Iron Endmembers from Aviris:

Red: Hematite	Green: Goethite	Blue: Jarosite
(high pH)		(low pH)

Note on Partial Unmixing:

”Unmixing is critical in imaging spectrometry since virtually every pixel contains some macroscopic mixture of materials. The theory of, and methods for unmixing of spectroscopic signatures are found in a number of sources. Hapke (1993) provides models for linear and non-linear spectral mixing and a discussion of the criteria for using each approach. In Earth remote sensing, a linear mixing model is typically used. Full unmixing which assumes that the spectra of all possible components are known, is described by van der Meer (2000) and Boardman (1989). Since it is often impossible to identify, a priori, all possible components, partial unmixing is an important tool. Match Filtering (MF) (Harsanyi and Chang 1994), Constrained Energy Minimization (CEM – similar to Match Filtering) (Farrand and Harsanyi 1997; Farrand 2001), and Mixture Tuned Match Filtering (${\text{ MTMF}}^{\text{ tm}}$) (Boardman et al. 1995) are commonly used methods for partial unmixing. ${\text{ MTMF}}^{\text{ tm}}$is probably the most popular unmixing method used for geologic remote sensing.”[2, Coulter sec 1.4.4]

(Emphasis Added.)

2 Cross-Correlation Spectral Matching

Cross-Correlation Spectral Matching

• $\{x_i;i=1...N\}$ wavelength vector with $N$ bands
• $\{y_i=y(x_i);i=1...N\}$ unknown (observed) image spectrum at given pixel
• $\{z_i=z(x_i);i=1...N\}$ known reference spectrum;
  (to be matched against y)

Linear Correlation Coefficient (Pearson’s $r$):

See [4, eq. (13.7.1)]. Generalize: cross-correlation at position $m<N$:

Ref. [8, eq. (33)]. Measures of fit:

1.

$r_0$ relative to 1 is not particularly good, but

$$t=r\sqrt{\frac{N-2}{1-r^2}}$$

(3)

(see [4, eq. (13.7.5) ff.]) is distributed in null case as Student’s $A\left(t\right|N-2)$.

2.

${\chi }^2$ statistic:

$${\chi }^2=(1-r^2)\sum _{i=1}^N{(y_i-\overline{y})}^2$$

(4)

See [4, eq. (14.2.13)]. For non-uniform weights, weight the sums in eq. (1) by $1∕{\sigma }_i^2$. Then goodness-of-fit

$$Q=Q\left(\frac{N-2}{2},\frac{{\chi }^2}{2}\right)$$

(5)

([4, eq. (14.2.12)]) where the incomplete gamma function $Q(a,x)$ is

$$\begin{array}{rcll}Q(a,x)& =& \frac{\Gamma (a,x)}{\Gamma \left(a\right)}=\frac{1}{\Gamma \left(a\right)}{\int \; <!--nolimits-->}_x^{\infty }{e}^{-t}{t}^{a-1}dt& \text{(6)}\\ Q(a,0)& =& 1\text{ and }Q(a,\infty )=0& \text{(7)}\end{array}$$

See [4, eq. (6.2.3)]. If e.g. $Q>0.1$, then the goodness-of-fit is believable. If $Q$ is larger than, say, 0.001, the fit may be acceptable if the the errors are non-normal or have been moderately underestimated. If $Q<0.001$ then question the model and/or estimation proceedure eq. (2). If latter, consider robust estimation.

Other useful measures Skewness

1.

Skewness: cross-correlogram of perfect match is parabola:

$$Skew\left(r_m\right)=\frac{1}{M}\sum _{m=1}^M{(r_m-r_{-m})}^{\left(2\right)}\ge 0$$

(8)

2.

RMS Skewness error [8, eq. (34)]:

$$RMS=\sqrt{\frac{\sum _{m=0}^M{(r_m-r_m^{\prime })}^2}{M}}$$

(9)

where $r_m^{\prime }$ is auto-correlogram of z with itself, $M$ is number of match positions, $m$ is match number.

1.

Another possibly useful statistical measure, if $(y,z)$ is binormal and $N$ is moderately large $(N\ge 10)$, is Fisher’s z-transform $z=(1∕2)\ln <!--FUNCTION\; APPLICATION-->\left[\right(1+r)∕(1-r\left)\right]$. Then if you measure the correlation coefficients of the unknown signal y against two reference signals $(z_1,z_2)$, the significance of a difference between the two measured correlation coefficients $r_1$ and $r_2$ is given by the complementary error function [4, eq. (13.7.10)]

$$\text{ erfc}\left(\frac{|r_1-r_2|}{\sqrt{\frac{2}{N_1-3}+\frac{2}{N_2-3}}}\right)$$

(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 $(z_1,z_2)$.

2.

Non-parametric or Rank Correlation (Spearman’s Rank-Order Correlation Coefficient). See [4, Sec. 13.8] and [6, ].

3 Optimized Cross-Correlation Spectral Matching

Optimized Cross-Correlation Spectral Matching In ”traditional” CCSM, the reference spectrum z to which $y$ is compared is taken to be a single (pure) endmember.
Next assume M endmembers $X_k$ of interest, and linear mixing. Seek weighting factors $\{a_k;k=1,\dots <!--FUNCTION\; APPLICATION-->M\}$

(synthesized pixel intensity at band $i$) that maximize the cross correlogram

at match position $m=0$, places that maximum at $m=0$ , and minimizes its skew, subject to $0\le a_k\le 1$ and ${\sum <!--FUNCTION\; APPLICATION-->}_{k=1}^M{a}_k=1$. Significance of different values of $r_0$ resulting from different choices of endmember set $\{X_k;k=1,\dots <!--FUNCTION\; APPLICATION-->M\}$ may be assessed using eq. (10). Coulter [2, ] chooses endmember set that maximizes $r_0$. As written (12) is the cross-correlation of two unit vectors:

which is independent of the normalization of $a$. Likewise, if $\overrightarrow{y}$ is spectrum at a particular image pixel, then the unit normalization $\widehat{y}$ rends the $r_m$ somewhat independent of shadow and striping.

Equation (14) makes $r_m$ non-linear in the $a_k$; eq (13) may be maximized w.r.t the $a_k$ 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 $y$. As 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 $r_m$ is expressed in terms of unit vectors, maximizing $r_0=\widehat{y}\cdot \widehat{z}$ is equivalent to minimizing

Minimizing the squared-difference of unit vectors $\widehat{y}$ and $\widehat{z}$ in (16) is still a non-linear problem in the $a_k$ because of how they appear in the denominator of the normalization of $\widehat{z}$ (14). However, if we relax the restrictions that z be a unit vector (and that $\sum <!--FUNCTION\; APPLICATION-->a_k\equiv 1$), we can define the equivalent problem

which is linear in the $a_k$, but not yet quite as general as we want. At each band $i$ there are (at least) three uncertainties in the sensor measurement: those in the received radiance $y_i$, those in the center wavelength $x_i$, 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 $x_i$ and assume normal independent distributions, we can write

as joint probability sensor records band radiance $y_i^0$ at recorded wavelength $x_i^0$, when it actually received (unknown) radiance $y_i$ at actual (and unknown) wavelength $x_i$. Joint probability of pixel’s measured spectrum across N bands is

Define

Associate the (unknown) actual spectrum $y\left(x_i\right)$ with the modeled mixture $z(x_i;a_1,\dots <!--FUNCTION\; APPLICATION-->a_M)={\sum <!--FUNCTION\; APPLICATION-->}_{k=1}^M{a}_k{X}_k\left(x_i\right)$. Then

represents a constrained optimization problem wherein we wish to minimize ${\chi }^2$ as a function of $N+M$ variables $(a_1,\dots <!--FUNCTION\; APPLICATION-->a_M,x_1,\dots <!--FUNCTION\; APPLICATION-->x_N)$ subject to $a_i\ge 0\forall <!--FUNCTION\; APPLICATION-->i$. The motivation for including the $(x_1,\dots <!--FUNCTION\; APPLICATION-->x_N)$ is to (hopefully) obtain ”more” non-negative $a_k$ in the optimal solution.

Normal Equations

where

and

and

Eqs (28) and (30) make eq. (27) quadratic in the $a_k$, so the optimization including center wavelengths uncertainties is non-linear. The $\frac{\partial X_j\left(x\right)}{\partial x}$ can be computed efficiently enough, but the $X_k\left(x_i\right)$ will need to be reconvolved against the new center wavelengths for each new set of center wavelengths $\{x_i;i=1\dots <!--FUNCTION\; APPLICATION-->N\}$ 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 $0$, equations (26) then become

which are the usual linear least-squares normal equations for the $a_k$. For what its worth, for the unconstrained problem the minimal solution vector $\{a_k;k=1,\dots <!--FUNCTION\; APPLICATION-->M\}$ will be unique provided the matrix $M_{ij}=X_i\cdot X_j$ is non-singular.

5 Image Classification and Statistics

5.1 Image Mean Vector and Covariance Matrix

Image Mean Vector and Covariance Matrix [9, White 2005] If image has $P$ bands and $N$ pixels, the mean vector is

where ${\overrightarrow{f}}_j$ is the jth pixel vector of the image

The image covariance matrix $C_{\text{ PxP}}$ is[1, Chang 2013 (6.6)]

where $F_{\text{ PxN}}$ is the matrix of $N$ pixel vectors each of length $P$

$M_{\text{ PxN}}$ is the matrix of $N$ identical mean vectors ($P$ rows by $N$ columns):

where $1_{\text{ PxN}}$ is an $P\times N$ matrix of ones.

5.2 Principal Component Transformation

Principal Component Transformation [5, Smith 1985] Karhunen-Loeve Transformation[10, White 2005] GRASS imagery i.pca: Let

$F_{\text{ PxN}}$

input-image multi-pixel vector ($P$ bands by $N$ pixels)

$M_{\text{ PxN}}$

mean vector matrix,

$Z_{\text{ PxN}}$

output-image multi-pixel vector,

$A_{\text{ PxP}}$

$P\times P$ matrix whose rows are the eigenvectors of the covariance matrix $C_{\text{ PxP}}$, arranged by decreasing magnitude of eigenvalue, as typically returned by SVD routines.

Magnitudes of ${\lambda }_i$ impose ordering on transformed component vectors $Z_{ik}={\sum <!--FUNCTION\; APPLICATION-->}_{j=1}^P{a}_{ij}{X}_{jk}$. Those with the largest ${\lambda }_i$ s.t. ${\lambda }_i∕{\lambda }_{max}>tol$ are the Principal Components. Tolerance tol should be related to the noise floor.

5.3 Minimum Noise Fraction

Minimum Noise Fraction[8, pg. 38] [3, ] We wish to find a particular coefficient matrix $\{a_{ij};i,j=1,\dots <!--FUNCTION\; APPLICATION-->P\}$ that in some sense maximizes the image S/N, assuming the image pixel vectors $\{{\overrightarrow{X}}_k,k=1,\dots <!--FUNCTION\; APPLICATION-->N\}$ are the sum of uncorrelated signal and noise:

Maximize

where ${\lambda }_i$ is generalized eigenvalue of $C_{\text{ PxP}}$ wrt $C_{\text{ PxP}}^N$, and $a_i$ are corresponding generalized eigenvectors. Compare with PCA:

Noise Covariance $C_{\text{ PxP}}^N$ Green suggests $C_{\text{ PxP}}$ be of unit variance and band-to-band uncorrelated.

• Sensor error estimate: (Aviris .rcc file) Call supplied error vector $\{e_j;j=1,\dots <!--FUNCTION\; APPLICATION-->P\}$. Then

  $$C_{lm}^N=e_m^n{\delta }_{lm}(n=1\text{ or }2)$$

  (56)

  is completely uncorrelated. In ideal case all $e_m$ are equal:

  $$\begin{array}{rcll}{C}_{lm}^N& =& eI,\text{ and }& \text{(57)}\\ a_i^T{C}_{\text{ PxP}}{a}_i^T& =& {\lambda }_i{a}_i^T{C}_{\text{ PxP}}^N{a}_i={\lambda }_{i}e\left(a_i^{T}I{a}_i\right)& \text{(58)}\\ & =& {\lambda }_{i}e(a_i^T\cdot a_i)={\lambda }_{i}e={\lambda }_i^{PCA}& \text{(59)}\end{array}$$

  since the eigenvectors $a_i$ are orthonormal. ${\lambda }_i={\lambda }_i^{PCA}$ if the variance $e\equiv 1$.

Homogeneous Area Method [7, sec 2.9.1] If possible find homogenous area of $N_h$ pixels in image:

Local Means and Local Variances [7, sec 2.9.2]

1.

Divide image into small $N_b$ pixel blocks (4x4, 5x5,...)

2.

For each (block, band) get local mean and variance: $$\begin{array}{rcll}{M}_{\text{ L}i}& =& \frac{1}{N_b}\sum _{k=1}^{N_b}{X}_{ik}i=1,2,\dots P\text{ bands}& \text{(65)}\\ {{\sigma }_{\text{ L}i}}^2& =& \frac{1}{N_b-1}\sum _{k=1}^{N_b}{(X_{ik}-M_{\text{ L}i})}^2& \text{(66)}\\ C_{\text{ L}ij}& =& \frac{\left({\delta }_{ij}\right)}{N_k-1}\sum _{k=1}^{N_b}(X_{ik}-M_{\text{ L}i})(X_{jk}-M_{\text{ L}j})& \text{(67)}\end{array}$$

$C$ is the local $N_b\times N_b$covariance matrix.

3.

bin $\left\{{{\sigma }_{\text{ L}i}}^2\right\}$ into classes between band min and max values.

4.

bin with most blocks represents mean noise of image. Hope this bin is same for all bands...

Local Means and Local Variances (con’t)

5.

Suppose ”most popular” bin is a P-cube, each side ranging $[{\sigma }_i^{\ast },{\sigma }_i^{\ast }+△{\sigma }_i^{\ast }]$, contains $N^{\ast }$ points. Then the average value over the bin

$$\overline{C_{\text{ L}ij}}=\frac{1}{N^{\ast }}\sum _{k=1}^{N^{\ast }}{\left(C_{\text{ L}ij}\right)}_k$$

(68)

is desired noise covariance matrix.

6.

Caveat: Assumes image is ”slowly varying enough” that enough of the $N_{\text{ total}}$ 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 small!

Other methods: ”unsupervised training” derived-endmember classification schemes e.g. LAS’ search ([11, ]) and GRASS’ cluster/maxlik are based upon local covariance minimization.

References