Application of the Stochastic Mixing Model to ... - Semantic Scholar
1924
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
Application of the Stochastic Mixing Model to Hyperspectral Resolution Enhancement Michael T. Eismann, Member, IEEE, and Russell C. Hardie, Senior Member, IEEE
Abstract—A maximum a posteriori (MAP) estimation method is described for enhancing the spatial resolution of a hyperspectral image using a higher resolution coincident panchromatic image. The approach makes use of a stochastic mixing model (SMM) of the underlying spectral scene content to develop a cost function that simultaneously optimizes the estimated hyperspectral scene relative to the observed hyperspectral and panchromatic imagery, as well as the local statistics of the spectral mixing model. The incorporation of the stochastic mixing model is found to be the key ingredient for reconstructing subpixel spectral information in that it provides the necessary constraints that lead to a well-conditioned linear system of equations for the high-resolution hyperspectral image estimate. Here, the mathematical formulation of the proposed MAP method is described. Also, enhancement results using various hyperspectral image datasets are provided. In general, it is found that the MAP/SMM method is able to reconstruct subpixel information in several principal components of the high-resolution hyperspectral image estimate, while the enhancement for conventional methods, like those based on least squares estimation, is limited primarily to the first principal component (i.e., the intensity component). Index Terms—Hyperspectral, Maximum a posteriori (MAP) estimation, resolution enhancement, stochastic mixing model.
I. INTRODUCTION
H
YPERSPECTRAL imaging is emerging as a valuable remote sensing technology for a variety of civilian, commercial, and military remote sensing applications. By exploiting fine spectral differences between various natural and man-made materials of interest, it can support improved detection and classification capabilities relative to panchromatic and multispectral remote sensors. A disadvantage of hyperspectral sensing, however, is that the spatial resolution is often coarser than panchromatic sensors. This engineering trade-off between spectral and spatial resolution arises in the sensor design as a result of the need to maintain imaging sensitivity at the finer spectral resolution of the hyperspectral sensor. Therefore, the enhanced spectral fidelity often comes at the expense of spatial fidelity, and fine shape and textural features can be lost. One approach that sensor designers have taken to address this issue is to couple low spatial resolution hyperspectral
Manuscript received October 10, 2003; revised April 10, 2004. M. T. Eismann is with the Sensors Directorate, Air Force Research Laboratory, Wright-Patterson AFB, OH 45433-7700 USA (e-mail: [email protected]). R. C. Hardie is with the Department of Electrical and Computer Engineering and the Electro-Optics Program, University of Dayton, Dayton, OH 45459-0226 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2004.830644
sensors with higher spatial resolution panchromatic and/or multispectral sensors. For example, the National Aeronautics and Space Administration (NASA) Earth Observing 1 satellite includes a 30-m hyperspectral sensor (Hyperion) and a 10-m panchromatic imager (Advanced Land Imager), and other commercial panchromatic imagers with resolution approaching 1 m are operational (e.g., Ikonos). Airborne hyperspectral sensors such as the Night Vision Imaging Spectrometer (NVIS) are also often flown with boresighted higher resolution panchromatic sensors. These sensor system approaches provide an opportunity to jointly process the hyperspectral and high-resolution imagery to achieve improved detection and/or classification performance. Joint processing of multiresolution imagery has been previously researched [1]–[9] in the context of improving the spatial detail in multispectral imagery, such as that from the Landsat Thematic Mapper, when coupled with higher resolution panchromatic imagery. The techniques that have been developed, however, are optimized to enhance visual information for human interpretation and are not particularly well suited to support automated processing associated with hyperspectral imagery. The primary difficulty is that the resolution enhancement problem is fundamentally ill-conditioned. Suppose that pixels and the low-resolution hyperspectral image contains spectral bands, and the high-resolution panchromatic image pixels where , usually by more than an contains pixels order of magnitude. The objective is to estimate the spectral bands of the high-resolution hyperspectral image ( unknowns) from the observed data ( linear , (256 equations). For example, let 256), and (64 64). Then, there are 14 680 064 unknowns to estimate from 983 040 equations. The ill-conditioned nature of the problem is usually addressed by constraining the available solutions. Principal component substitution (PCS) does this by only allowing the first principal component of the estimate to contain the high-resolution content of the panchromatic image. Since the first principal component of hyperspectral data generally corresponds to the intensity component, this only allows the intensity of the resulting hyperspectral image to vary at the finer resolution level, and not the spectral characteristics. Least squares estimation (LSE), a methodology upon which some of the prior resolution enhancement methods are based or to which they are related, constrains the estimated high-resolution image spectra corresponding to each low-resolution image spectrum to a line in the multidimensional space determined by the local correlation characteristics between the low-resolution hyperspectral image and the panchromatic image. The resulting resolution enhancement of
0196-2892/04$20.00 © 2004 IEEE
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
1925
this approach is also generally limited to the first principal component. This is illustrated later in this paper, and a more intuitive description is given in Section II. Recent research in hyperspectral resolution enhancement has focused on the incorporation of a spectral mixture model into the estimation problem. One approach starts with conventional linear unmixing to generate abundance maps, and then uses constrained optimization based on the panchromatic image and a strict radiometric model to spatially locate the endmembers to high resolution [10]–[12]. Published results are limited to well-behaved synthetic data where the endmembers are known a priori. An alternative approach appends the high-resolution image to the hyperspectral data and computes a mixture model based on the joint dataset [13]. A high-resolution hyperspectral image, however, is not explicitly estimated. In this paper, we describe the application of a maximum a posteriori (MAP) estimation framework [14] with a stochastic mixing model (SMM) to the hyperspectral resolution enhancement problem. The stochastic mixing model provides a physical model to constrain the estimation problem in a similar manner as the linear mixing model in the approaches described in [10]–[12]. However, by using a statistical, as opposed to deterministic, framework, the approach does not rely on strict radiometric conformance and should be less sensitive to sensor degradations (noise, miscalibration, etc.) that reduce the radiometric correspondence between the two observations. The remainder of the paper is organized as follows. In Section II, the mathematical formulation of the MAP/SMM estimation approach is described, including comparison to the PCS and LSE methods (which will serve as performance benchmarks). In Section III, the characteristics of the test imagery and the quantitative metrics used to assess resolution enhancement performance are given. Experimental results are presented in Section IV, and conclusions are given in Section V.
of hyperspectral image pixels. This is a band-interleaved-bypixel spatially rastered representation of the inherently threedimensional (3-D) data. The objective of hyperspectral resolution enhancement is to estimate (or reconstruct) what is termed the high-resolution hyperspectral image, i.e., one possessing the spectral content of the low-resolution hyperspectral image but at the spatial resolution of the panchromatic image. This image is denoted as
II. MATHEMATICAL FORMULATION Suppose that an area is remotely imaged concurrently with a hyperspectral imaging system and a higher resolution panchromatic imaging system. The hyperspectral image observation will be termed “low-resolution” from the standpoint that its spatial resolution is lower than the panchromatic image. The panchromatic image is described as a one-dimensional (1-D) vector of intensity measurements (1) is the panchromatic image intensity at the spatial where is the number of location designated by the index , and panchromatic image pixels. This 1-D, or lexicographic, format is a rastered representation of the two-dimensional image. In a similar manner, the low-resolution hyperspectral image is also represented as a 1-D vector of spectroradiometric measurements
(3) where is a high-resolution pixel spectrum (with bands) at the spatial location designated by index , and is the number of image pixels (both corresponding to the panchromatic image). This is again a band-interleaved-by-pixel spatially rastered representation of the inherently 3-D data. A. Observation Model The relationship between the panchromatic image and the high-resolution hyperspectral image can be written in matrix form as (4) where the spectral response matrix is a sparse matrix whose columns are the spectral response functions for the panchromatic pixel locations, and is assumed to be a spatially independent zero-mean Gaussian random process with a standard deviation of . Similarly, the low-resolution hyperspectral image is related to the high-resolution hyperspectral image by (5) where is a sparse matrix whose rows are the spatial response functions for the low-resolution hyperspectral pixels and is assumed to be a spatially independent zero-mean Gaussian . Note random process with a spectral covariance matrix that both and are noninvertible, nonsquare matrices, so there is no unique solution for in either equation even in the absence of observation noise. B. MAP Formulation The goal of MAP estimation, as applied to the hyperspectral resolution problem, is to find an estimate for the high-resolution hyperspectral image that maximizes its conditional probability relative to the two observations, or (6) Using Bayes’ rule, the conditional probability density function can be expressed as (7)
(2) where is the low-resolution pixel spectrum (with bands) at the spatial location designated by index , and is the number
If the random noise processes of the two observations are not correlated (generally a good assumption given that the panchromatic and hyperspectral data come from different sensors), then
1926
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
the observed images are conditionally independent and can be written as the product of two conditional density functions (8)
Note that this function is quadratic in , so that a closed-form solution for the minimum can be readily obtained by determining the where the gradient of the cost function is zero, and its second derivative is positive. With some mathematical manipulation, it can be shown that this estimate corresponds to the solution to the linear system of equations
Ignoring the term in the denominator, which does not depend on , the MAP estimate is then given by (16) (9) The conditional density functions are easily derived from the observation models discussed in Section I. They take the form of Gaussian density functions by virtue of the assumed noise statistics
where (17) and (18)
(10) and
The Hessian matrix is a square, positive semidefinite matrix. Therefore, a solution to (16) will always represent a minimum. D. Stochastic Mixing Model
(11) The density function for the high-resolution hyperspectral image is unknown. Assume, however, that it is normally distributed. Then, it can be expressed as
(12) where the mean vector and covariance matrix can, in general, vary as a function of spatial location. Assuming spatial indepenwill be block diagonal. More specifically, the mean dence, vector and covariance matrix will take the forms (13) and
The solution of (16) requires some underlying model for the high-resolution hyperspectral image, or more specifically, the mean vector and covariance matrix as a function of spatial position. The stochastic mixing model is derived from the low-resolution hyperspectral image and used for this purpose. The formulation and implementation of the stochastic mixing model is described in this subsection. For a more detailed description, see the basic formulation in [16] and refinements in [17]. The stochastic mixing model is similar to the well-developed linear mixing model in that it attempts to decompose spectral data in terms of a linear combination of endmember spectra. However, in the case of the stochastic mixing model, the endare normally distributed random vectors, parammembers and covariance matrix eterized by their mean vector , instead of deterministic spectra. The endmember index ranges from one to , the number of endmembers assumed of the observed hyperspectral in the scene. Each spectrum image is then assumed to be a member of a mixture class defined by the random vector , where
(14) where
(19)
refers to a direct sum [15].
C. Optimization From the form of the probability density functions that compose (9), an objective function can be formed whose minimization will result in the optimal estimate. This objective function is given by the exponents of (10)–(12)
(15)
are the endmember abundances associated with the mixture class, and is the mixture class index. The endmember abundances are assumed to conform to the physical constraints of being between zero and one, and summing to unity, and are quantized at a specified level. This is done by a structured, recursive search algorithm that determines all the possible combinations of quantized abundances (i.e., discrete levels between zero and one) for the specified number of endmembers that satisfy the physical constraints. The result is a finite and discrete set of mixture classes.
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
1927
Due to the linear relationship in (19), the mean vector and covariance matrix for the th mixture class are functionally dependent on the corresponding endmember statistics by (20) and (21) and the class-conditional probability density function is described by the Gaussian mixture distribution Fig. 1. Graphical depiction of conventional LSE method.
(22) where
(23) The mixture model parameters that must be estimated from the measured data are the prior probabilities for all the mixture classes and the mean vectors and covariance matrices of the endmember classes. The discrete set of feasible classes are prede. fined by the set of abundances The stochastic expectation maximization [18] approach is used to self-consistently estimate the model parameters , , and , and the mixture class assignments [i.e., abundance maps ] for the observed low-resolution hyperspectral scene data. The specific implementation used in this paper incorporates the modifications detailed in [17]. It proceeds in the following manner. First, the endmember statistics are initialized. The mean vectors are initialized to a set of well-separated image spectra selected from the image by simplex maximization [19], the covariance matrices are initialized to the global covariance matrix, and the prior probabilities are uniformly distributed. Then, a series of iterations (nominally 20) are performed, each consisting of the following three steps. Step 1) The posterior probability is estimated on the th iteration for all combinations of image spectra and mixture classes according to (24)
Step 2) The image spectra are classified into the mixture classes in (19) by a Monte Carlo process based on the posterior probability values. Step 3) The endmember mean and covariance estimates are updated based on the sample statistics of spectra classified into the respective pure classes, and the prior probability estimates are updated based on the relative populations of the mixture classes.
After the algorithm converges, the low-resolution image spectra are ultimately assigned to the mixture classes , , that maximize the conditional posterior probability in (24), and low-resolution abundance maps are formed. That is, there are abundance maps giving the spatial abundance distributions for each endmember. The estimated abundance maps are then bilinearly , interpolated to give the high-resolution abundance maps , . According to (20) and (21), the elements of the spatially varying statistics in (13) and (14) are then estimated by (25) and (26) This provides the remaining terms in (17) and (18) to solve the system of equations in (16) for the high-resolution hyperspectral image estimate. E. Comparison to LSE In this paper, an LSE solution to the resolution enhancement problem is used as a performance benchmark. The LSE estimate takes the form (27) is the average of those panchromatic intensity values where that span the portion of the scene corresponding to . The vector can be estimated as a cross covariance [6] or using a least squares linear regression technique [1] based on and the available lower resolution images (e.g., the set of samples). The cross-covariance approach is corresponding used in this paper. A graphical comparison of the LSE and MAP/SMM methods is given in Figs. 1 and 2 for a simple case consisting of a two and corresponding to a panchromatic intensity values single, two-band, low-resolution spectrum . The LSE method constrains the estimated high-resolution spectra and to the
1928
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
III. DATA AND METRICS A. Test Data
Fig. 2.
Graphical depiction of proposed MAP/SMM method.
indicated line given by (27). The estimates are symmetrically located about (i.e., their average spectrum equals ) with a separation proportional to the difference in panchromatic intensity values. The MAP/SMM method statistically balances consistency with the observation models for and along with the mixture class probability density function from the interpolated abundance maps and SMM endmember statistics. This is depicted in Fig. 2 for the same simple case and uniform spatial and spectral response functions. Consistency with the model for forces the estimated high-resolution spectra toward the two parallel lines in the figure while consistency with occurs when the estimates are symmetrically located about . The SMM forces the estimates into their respective mixture classes. Variation can result toward the corresponding endmember classes and occurs particularly when neighboring pixels are purer (due to the linear interpolation of the abundance maps). F. Computational Issues The Hessian matrix is an matrix, so (14) represents a very large system of equations. It is also not necessarily a well-conditioned system of equations because of the high degree of correlation that often exists in hyperspectral data. That is, the rank of the covariance matrices associated with hyperspectral imagery is often less than the number of spectral bands. This issue is addressed by processing in the principal components domain rather than the spectral domain and limiting the data to a leading order subspace corresponding to the highest eigenvalues. Most hyperspectral imagery can be well represented with 10–20 principal components. In this paper, we also make the assumption that the point spread function corresponding to neighboring low-resolution pixels is nonoverlapping. Optically, this corresponds to the case where the resolution of the low-resolution hyperspectral image is strictly defined by system of equations the detector area. It allows the of smaller ( to be decomposed into a number ) systems equations, which significantly reduces the required computation. Computationally efficient methods to solve the problem without employing this assumption, however, are currently being studied.
The resolution enhancement results presented in this paper are based on three hyperspectral data sources. The first is a synthetically generated image, which provides a well-controlled test case where the imagery is known to conform to the basic model assumptions. The image was generated by the Digital Imaging and Remote Sensing Synthetic Image Generator (DIRSIG) model [20]. DIRSIG is a rigorous end-to-end spectroradiometric scene model that incorporates full 3-D scene construction; radiation propagation, absorption, and scattering; and extensive target-background interactions. The particular scene used is an emulation of data from the Hyperspectral Digital Imagery Collection Experiment (HYDICE) airborne sensor [21] viewing a forested area containing a clearing with a dirt road. The spatial resolution is approximately 0.15 m, about four to six times better than that collected by the HYDICE sensor during low-altitude operations (1.5–2-km altitude with 0.5-mrad angular resolution). The spectral range (0.38–2.6 m) and resolution (ranging from 4–20 nm) are based on the HYDICE sensor characteristics. The second image was collected by the Airborne Visible-Infrared Imaging Spectrometer (AVIRIS) sensor [22]. AVIRIS is a scanning dispersive hyperspectral imaging sensor that flies on the NASA ER-2 aircraft nominally 20 km above ground level and covers the visible/near-infrared (VNIR) and shortwave-infrared (SWIR) spectral regions. The spatial resolution is about 20 m, the spectral range is from 0.37–2.5 m, and the spectral resolution is on the order of 10 nm. The image used in this study was collected over Yorktown, VA. Because of the degree of complexity and low resolution of the image, it poses a challenge for resolution enhancement. The third image was collected by the NVIS airborne sensor [23]. NVIS is a pushbroom dispersive imaging sensor that flies on the Twin Otter aircraft nominally 1 km above ground level and also covers the VNIR/SWIR spectral range. The spatial resolution is about 0.6 m, the spectral range is from 0.45–2.5 m, and the spectral resolution is on the order of 5–10 nm. The image used in this study was collected over Fort A.P. Hill, VA, on May 4, 2000. The analysis presented in this paper is based only on the VNIR hyperspectral data covering the 0.45–0.9 m spectral range. In the DIRSIG and AVIRIS cases, the original imagery was spatially degraded by a factor of four to produce the low-resolution hyperspectral observation, and spectrally integrated over the entire spectral range to produce the high-resolution panchromatic observation. This produces imagery representative of hypothetical sensors with poorer spatial resolution and higher SNR (due to the spatial averaging) relative to the original imagery. This means that results obtained using this approach are not representative of data from the AVIRIS sensor were it to be used directly for the low-resolution hyperspectral data. Thus, the quantitative conclusions drawn are limited to comparing the MAP/SMM and LSE methods, and not the achievable enhancement for realistic AVIRIS data. The NVIS sensor also contains a bore-sighted VNIR panchromatic sensor of higher spatial resolution (0.5 ft), making it a good test set for the realistic employment of the resolution
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
Fig. 3. Panchromatic image for DIRSIG test case.
Fig. 4. Panchromatic image for AVIRIS test case.
1929
spatial and spectral characteristics, the metrics used to quantify the performance are designed to capture the difference in these fine details. Quantitative analysis of the NVIS test case is not possible because the lack of a truth image. The rms difference between the high-resolution hyperspectral image estimate and the truth image is the most basic error metric used. It is computed on a band-by-band basis. The spectral angle difference between the estimate and the truth for each spatial position is also computed. This metric isolates the errors in the spectral shape of the estimate from (correlated) intensity errors. In this paper, histograms of the spectral angle difference are used to compare methods. Neither of these metrics directly measures the ability of the hyperspectral resolution enhancement methods to reconstruct spatial content in the imagery beyond the resolution of the lowresolution hyperspectral sensor, the primary objective of the algorithms. To quantify this ability, a metric was developed that measures the correlation of the subpixel image structure to the truth. This is performed in the following manner. First, the estimate and the truth image are transformed to the principal component domain using the eigenvector matrix associated with the low-resolution spectral covariance matrix. Next, the mean intensity for each set of spatial pixels associated with a low-resolution pixel (i.e., super-pixel) is computed and subtracted. This is performed independently for each principal component band. Finally, the subpixel correlation metric is computed on a band-byband basis as the ratio of the cross correlation between the highpass filtered estimate with the truth image relative to the autocorrelation of the truth image
(28)
Fig. 5. Panchromatic image for NVIS test case.
enhancement method. Spatial registration of the panchromatic to hyperspectral imagery was performed by a three-step process consisting of local phase correlation measurement of relative image misregistration, spline interpolation of the measurements, and resampling with a finite-impulse-response interpolation kernel [24], [25]. The panchromatic imagery for the three datasets are illustrated in Figs. 3–5. B. Performance Metrics Four metrics were developed to compare the estimate produced by the resolution enhancement methods to the high-resolution truth image. Because the difference between the resolution enhancement methods lies in subtle, but important, fine
is the filtered estimate for the th spatial position and where the th principal component, and is the filtered truth for the same position and principal component. In the case of the DIRSIG forest scene, a spectral matched filter was constructed to detect the dirt road, a fine-resolution image feature that serves as a fourth test of the resolution enhancement methods. The filter vector was computed as the mean spectrum over the first 50 pixels in column 164, aligned directly on a soil track in the dirt road. Because this image is dominated by vegetation, which has fairly homogeneous spectral content, the spectral matched filter works very well at detecting the soil feature in the high-resolution image. The separation of the two soil tracks, however, is not resolved in the low-resolution result. Receiver operating characteristic (ROC) curves were generated in the following manner. First, a truth mask was made by thresholding the spectral matched filter output for the true high-resolution image at a level that visibly captured the soil track feature. A matched filter output was generated for each high-resolution image estimate and thresholded at a range of levels. The outputs above threshold at each level were then compared to the truth mask and scored in terms of the fraction detected and false-alarm rate to form the ROC curve.
1930
Fig. 6.
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
RMS difference for DIRSIG test case.
Fig. 7. Spectral angle error for DIRSIG test case.
IV. RESULTS A. DIRSIG Imagery Resolution-enhanced hyperspectral imagery using the PCS, LSE, and MAP/SMM methods was produced and quantitatively analyzed using the aforementioned metrics. These results are given in Figs. 6–9, where the SMM consisted of four endmem) and 74 mixture classes (i.e., ). There bers (i.e., is a modest improvement in the rms error relative to LSE with the MAP/SMM method, and a substantial improvement relative to PCS. Note that wavelet-based merging of the panchromatic image into the first principal component [26] is likely to produce a result similar to LSE. Similar trends occur for the spectral angle error results, where the MAP/SMM method provides a significant reduction in the number of pixels with spectral angle error in the 10 to 30 range relative to LSE. The subpixel cross-correlation results support the hypothesis that incorporation of the SMM will support the accurate reconstruction of subpixel spectral content, even down to the lower principal components. Enhancement for the LSE method is limited to the first
Fig. 8. Subpixel cross correlation for DIRSIG test case.
Fig. 9.
ROC curve for soil track detectability for DIRSIG test case.
principal component. Because detection of the soil track feature relies on high-resolution spectral information, the ROC curve associated with this detection problem is significantly improved by the MAP/SMM resolution enhancement. This is probably the best indicator of the efficacy of the MAP/SMM approach for this dataset. Qualitatively, the MAP/SMM method was found to be able to reconstruct subpixel spatial structure in the lower order principal components, while the LSE method is only able to enhance the first principal component. There was no visible difference between the first principal component imagery for both cases, but there were substantial differences in the lower order components. For example, Fig. 10 compares the second principal components for both cases as compared to the low-resolution observation and high-resolution truth. B. AVIRIS Imagery The results based on the AVIRIS imagery were fairly consistent with the DIRSIG results. Again, the low-resolution SMM
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
(a)
1931
(b)
Fig. 12. Subpixel cross correlation for AVIRIS test case.
(c)
(d)
Fig. 10. Second principal component images for DIRSIG test case. (a) Low resolution. (b) LSE estimate. (c) MAP/SMM estimate. (d) High-resolution truth.
Fig. 11.
Spectral angle error for AVIRIS test case.
(a)
(b)
(c)
(d)
Fig. 13. Third principal component images for AVIRIS test case. (a) Low resolution. (b) LSE estimate. (c) MAP/SMM estimate. (d) High-resolution truth.
consisted of four endmembers and 74 mixture classes. A subset of the quantitative and qualitative results is given in Figs. 11–13. Similarly to the DIRSIG case, the MAP/SMM method was able to reconstruct fine-resolution structure in the lower order principal components that LSE was not. This is clearly evident in the magnified, third principal component images in Fig. 13. Relative to the DIRSIG results, the quantitative improvement in spectral angle error and subpixel cross correlation compared to LSE is more modest for the AVIRIS test case.
only the first principal component contained significantly enhanced spatial resolution relative to the low-resolution image (not shown). The lower principal components were visibly identical to the low-resolution image. There is, however, higher resolution content in all principal components images for the MAP/SMM result. In this case, there was no high-resolution image truth to provide a quantitative comparison. Instead, the purpose was to demonstrate the efficacy of the algorithm in the real situation, i.e., coincident imagery from separate sensors.
C. NVIS Imagery
D. Sensitivity Analysis
Again, four endmembers and 74 mixture classes were used to derive the SMM for the NVIS case. The resulting principal component images are compared in Fig. 14. For the LSE result,
A sensitivity analysis was performed for both the DIRSIG and AVIRIS test cases to assess the sensitivity of MAP/SMM resolution enhancement performance to panchromatic image noise
1932
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
(a)
Fig. 15. Soil track detectability versus panchromatic image SNR for DIRSIG test case.
(b)
(c)
Fig. 16. Soil track detectability versus hyperspectral image resolution factor for DIRSIG test case.
(d) Fig. 14. Principal component images for NVIS test case. (Left) LSE and (right) MAP/SMM. (a) First principal component. (b) Second principal component. (c) Third principal component. (d) Fourth principal component.
and hyperspectral image spatial resolution. The results are presented here in terms of the DIRSIG test case and the soil track detectability metric (the most sensitive metric), but the trends for the AVIRIS imagery and the other metrics were consistent. Fig. 15 shows the detectability results for a range of panchromatic image SNRs, defined as the ratio of the image-to-noise standard deviation. There is no reduction in performance until
the SNR drops into the 10–100 range, but the detectability drops quickly below an SNR of 10. Fig. 16 illustrates the reduction in performance as the spatial resolution of the hyperspectral observation image is varied while the panchromatic image is not changed. The legend indicates the ratio of panchromatic to hyperspectral spatial resolutions. The baseline case is a resolution factor of 4. A steady decrease in detectability with hyperspectral image spatial resolution is observed. V. CONCLUSION The results presented in this paper have provided an experimental proof-of-concept that the incorporation of a stochastic mixing model derived from the low-resolution hyperspectral observation in a maximum a posteriori estimator can provide better hyperspectral resolution enhancement results than conventional methods such as principal component substitution and
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
1933
least squares estimation. For all of the datasets investigated, the MAP/SMM method was able to produce more accurate subpixel structure in the lower order principal component images than the conventional methods. Additionally, the method was found to be robust to variations (e.g., noise and spatial resolution) in the observation imagery. We anticipate that such resolution enhancement will have positive impacts on subsequent hyperspectral processing and analysis, especially where finer spatial resolution is desired. This, however, requires further investigation. Other areas where further work is warranted include generalizing the implementation to the case of overlapping point spread functions, incorporating multispectral high-resolution imagery, and thoroughly analyzing the noise characteristics of the enhanced estimate relative to the low-resolution observation. We expect some inherent reduction in SNR with resolution enhancement, but have not yet characterized this effect.
[14] R. C. Hardie and M. T. Eismann, “MAP estimation for cross-sensor spatial resolution enhancement,” IEEE Trans. Image Processing, vol. 13, pp. 1174–1184, Sept. 2004. [15] H. Eves, Elementary Matrix Theory. New York: Dover, 1966, p. 107. [16] A. D. Stocker and A. P. Schaum, “Application of stochastic mixing models to hyperspectral detection problems,” Proc. SPIE, vol. 3071, pp. 47–60, 1997. [17] M. T. Eismann and R. C. Hardie, “Initialization and convergence of the stochastic mixing model,” Proc. SPIE, vol. 5159, pp. 307–318, 2003. [18] P. Masson and W. Pieczynski, “SEM algorithm and unsupervised statistical segmentation of satellite images,” IEEE Trans. Geosci. Remote Sensing, vol. 31, pp. 618–633, May 1993. [19] M. E. Winter, “Fast autonomous endmember determination in hyperspectral data,” in Proc. 13th Int. Conf. Applied Geological Remote Sensing, vol. II, 1999, pp. 337–344. [20] J. R. Schott et al., “Incorporation of time-dependent thermodynamic model and radiation propagation model into infrared three-dimensional synthetic image generation,” Opt. Eng., vol. 37, no. 7, p. 1505, July 1992. [21] R. R. Basedow et al., “The HYDICE instrument design,” in Proc. ISSR, Nov. 1992. [22] R. O. Green et al., “Imaging spectroscopy and the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS),” Remote Sens. Environ., vol. 65, pp. 227–248, 1998. [23] C. Simi et al., “Night Vision Imaging Spectrometer (NVIS) performance parameters and their impact on various detection algorithms,” Proc. SPIE, vol. 4049, pp. 218–229, 2000. [24] A. Schaum and M. McHugh, “Analytic methods of image registration: Displacement estimation and resampling,” Nav. Res. Lab., West Bethesda, MD, NRL Rep. 9298, 1991. [25] P. Landcaster and K. Salkauskas, Curve and Surface Fitting: An Introduction. London, U.K.: Academic, 1986. [26] T. Ranchin and L. Wald, “Fusion of high spatial and spectral resolution images: the ARSIS concept and its implementation,” Photogram. Eng. Remote Sens., vol. 66, pp. 49–56, 2000.
ACKNOWLEDGMENT The authors would like to acknowledge the support of R. Bartell (Veridian Corporation), B. Kendall (Space Computer Corporation), E. Winter (Technical Research Associates), and the NASA Jet Propulsion Laboratory for the imagery used in this paper. REFERENCES [1] J. C. Price, “Combining panchromatic and multispectral imagery from dual resolution satellite instruments,” Remote Sens. Environ., vol. 21, pp. 119–128, 1987. [2] C. K. Munechika, J. S. Warnick, C. Salvaggio, and J. R. Schott, “Resolution enhancement of multispectral image data to improve classification accuracy,” Photogram. Eng. Remote Sens., vol. 59, pp. 67–72, 1993. [3] W. J. Carper, T. M. Lillesand, and R. W. Kiefer, “The use of intensity-hue-saturation transformations for merging SPOT panchromatic and multispectral image data,” Photogram. Eng. Remote Sens., vol. 56, pp. 459–467, 1990. [4] V. K. Shettigara, “A generalized component substitution technique for spatial enhancement of multispectral images using a higher resolution data set,” Photogram. Eng. Remote Sens., vol. 58, pp. 561–567, 1992. [5] A. E. Iverson and J. R. Lersch, “Adaptive image sharpening using multiresolution representations,” Proc. SPIE, vol. 2231, pp. 72–83, 1994. [6] R. Nishii, S. Kusanobu, and S. Tanaka, “Enhancement of low resolution image based on high resolution bands,” IEEE Trans. Geosci. Remote Sensing, vol. 34, pp. 1151–1158, Sept. 1996. [7] D. P. Filberti, S. E. Marsh, and R. A. Schowengardt, “Synthesis of imagery from high spatial and spectral resolution from multiple image sources,” Opt. Eng., vol. 33, no. 8, Aug. 1994. [8] B. Zhukov, M. Berger, F. Lanzl, and H. Kaufman, “A new technique for merging multispectral and panchromatic images revealing sub-pixel spectral variation,” in Proc. IGARSS, Firenze, Italy, July 10–14, 1995, pp. 2154–2156. [9] B. Zhukov, D. Oertel, and F. Lanzl, “A multiresolution, multisensor technique for satellite remote sensing,” in Proc. IGARSS, Firenze, Italy, July 10–14, 1995, pp. 51–53. [10] H. H. Gross and J. R. Schott, “Application of spectral mixture analysis and image fusion techniques for image sharpening,” Remote Sens. Environ., vol. 63, pp. 85–94, 1998. [11] G. D. Robinson, H. N. Gross, and J. R. Schott, “Evaluation of two applications of spectral mixing models to image fusion,” Remote Sens. Environ., vol. 71, pp. 272–281, 2000. [12] H. N. Gross and J. R. Schott, “Application of spatial resolution enhancement and spectral mixture analysis to hyperspectral images,” Proc. SPIE, vol. 2821, pp. 30–41, 1996. [13] M. E. Winter and E. M. Winter, “Physics-based resolution enhancement of hyperspectral data,” Proc. SPIE, vol. 4725, pp. 580–587, 2002.
Michael T. Eismann (M’03) received the B.S. degree in physics from Thomas More College, Crestview Hills, KY, in 1985, the M.S. degree in electrical engineering from the Georgia Institute of Technology, Atlanta, in 1987, and the Ph.D. degree in electrooptics from the University of Dayton, Dayton, OH, in 2004. He is the currently a Technical Advisor with the Electro-Optical Targeting Branch, Sensors Directorate, Air Force Research Laboratory (AFRL), WrightPatterson Air Force Base, OH. He is responsible for overseeing the development of passive electrooptical and infrared sensor technology, and transition into operational airborne targeting and reconnaissance systems. His recent research emphasis has been the development of hyperspectral imaging sensors and image processing. Prior to joining AFRL in 1996, he was with the Environmental Research Institute of Michigan (ERIM). He joined ERIM in 1991, where he was involved in research concerning active and passive optical and infrared targeting and reconnaissance, optical information processing, and holographic optics.
Russell C. Hardie (S’91–M’92–SM’99) received the B.S. degree (magna cum laude) in engineering science from Loyola College, Baltimore, MD, in 1988, and the M.S. and Ph.D. degrees in electrical engineering from the University of Delaware, Newark, in 1990 and 1992, respectively. He is currently an Associate Professor with the Department of Electrical and Computer Engineering, University of Dayton, OH, and holds a joint appointment with the Electro-Optics Program. He was a Senior Scientist with the Earth Satellite Corporation prior to his appointment at the University of Dayton in 1993. His research interests include a wide variety of topics in the area of digital signal and image processing, including image restoration and medical image processing. Dr. Hardie was a corecipient of the Rudolf Kingslake Medal and Prize from SPIE in 1998, for work on multiframe image resolution enhancement algorithms. In 1999, he received the School of Engineering Award of Excellence in Teaching at the University of Dayton. He received the First Annual Professor of the Year Award in 2002 from the student chapter of the IEEE at the University of Dayton.
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
Application of the Stochastic Mixing Model to Hyperspectral Resolution Enhancement Michael T. Eismann, Member, IEEE, and Russell C. Hardie, Senior Member, IEEE
Abstract—A maximum a posteriori (MAP) estimation method is described for enhancing the spatial resolution of a hyperspectral image using a higher resolution coincident panchromatic image. The approach makes use of a stochastic mixing model (SMM) of the underlying spectral scene content to develop a cost function that simultaneously optimizes the estimated hyperspectral scene relative to the observed hyperspectral and panchromatic imagery, as well as the local statistics of the spectral mixing model. The incorporation of the stochastic mixing model is found to be the key ingredient for reconstructing subpixel spectral information in that it provides the necessary constraints that lead to a well-conditioned linear system of equations for the high-resolution hyperspectral image estimate. Here, the mathematical formulation of the proposed MAP method is described. Also, enhancement results using various hyperspectral image datasets are provided. In general, it is found that the MAP/SMM method is able to reconstruct subpixel information in several principal components of the high-resolution hyperspectral image estimate, while the enhancement for conventional methods, like those based on least squares estimation, is limited primarily to the first principal component (i.e., the intensity component). Index Terms—Hyperspectral, Maximum a posteriori (MAP) estimation, resolution enhancement, stochastic mixing model.
I. INTRODUCTION
H
YPERSPECTRAL imaging is emerging as a valuable remote sensing technology for a variety of civilian, commercial, and military remote sensing applications. By exploiting fine spectral differences between various natural and man-made materials of interest, it can support improved detection and classification capabilities relative to panchromatic and multispectral remote sensors. A disadvantage of hyperspectral sensing, however, is that the spatial resolution is often coarser than panchromatic sensors. This engineering trade-off between spectral and spatial resolution arises in the sensor design as a result of the need to maintain imaging sensitivity at the finer spectral resolution of the hyperspectral sensor. Therefore, the enhanced spectral fidelity often comes at the expense of spatial fidelity, and fine shape and textural features can be lost. One approach that sensor designers have taken to address this issue is to couple low spatial resolution hyperspectral
Manuscript received October 10, 2003; revised April 10, 2004. M. T. Eismann is with the Sensors Directorate, Air Force Research Laboratory, Wright-Patterson AFB, OH 45433-7700 USA (e-mail: [email protected]). R. C. Hardie is with the Department of Electrical and Computer Engineering and the Electro-Optics Program, University of Dayton, Dayton, OH 45459-0226 USA (e-mail: [email protected]). Digital Object Identifier 10.1109/TGRS.2004.830644
sensors with higher spatial resolution panchromatic and/or multispectral sensors. For example, the National Aeronautics and Space Administration (NASA) Earth Observing 1 satellite includes a 30-m hyperspectral sensor (Hyperion) and a 10-m panchromatic imager (Advanced Land Imager), and other commercial panchromatic imagers with resolution approaching 1 m are operational (e.g., Ikonos). Airborne hyperspectral sensors such as the Night Vision Imaging Spectrometer (NVIS) are also often flown with boresighted higher resolution panchromatic sensors. These sensor system approaches provide an opportunity to jointly process the hyperspectral and high-resolution imagery to achieve improved detection and/or classification performance. Joint processing of multiresolution imagery has been previously researched [1]–[9] in the context of improving the spatial detail in multispectral imagery, such as that from the Landsat Thematic Mapper, when coupled with higher resolution panchromatic imagery. The techniques that have been developed, however, are optimized to enhance visual information for human interpretation and are not particularly well suited to support automated processing associated with hyperspectral imagery. The primary difficulty is that the resolution enhancement problem is fundamentally ill-conditioned. Suppose that pixels and the low-resolution hyperspectral image contains spectral bands, and the high-resolution panchromatic image pixels where , usually by more than an contains pixels order of magnitude. The objective is to estimate the spectral bands of the high-resolution hyperspectral image ( unknowns) from the observed data ( linear , (256 equations). For example, let 256), and (64 64). Then, there are 14 680 064 unknowns to estimate from 983 040 equations. The ill-conditioned nature of the problem is usually addressed by constraining the available solutions. Principal component substitution (PCS) does this by only allowing the first principal component of the estimate to contain the high-resolution content of the panchromatic image. Since the first principal component of hyperspectral data generally corresponds to the intensity component, this only allows the intensity of the resulting hyperspectral image to vary at the finer resolution level, and not the spectral characteristics. Least squares estimation (LSE), a methodology upon which some of the prior resolution enhancement methods are based or to which they are related, constrains the estimated high-resolution image spectra corresponding to each low-resolution image spectrum to a line in the multidimensional space determined by the local correlation characteristics between the low-resolution hyperspectral image and the panchromatic image. The resulting resolution enhancement of
0196-2892/04$20.00 © 2004 IEEE
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
1925
this approach is also generally limited to the first principal component. This is illustrated later in this paper, and a more intuitive description is given in Section II. Recent research in hyperspectral resolution enhancement has focused on the incorporation of a spectral mixture model into the estimation problem. One approach starts with conventional linear unmixing to generate abundance maps, and then uses constrained optimization based on the panchromatic image and a strict radiometric model to spatially locate the endmembers to high resolution [10]–[12]. Published results are limited to well-behaved synthetic data where the endmembers are known a priori. An alternative approach appends the high-resolution image to the hyperspectral data and computes a mixture model based on the joint dataset [13]. A high-resolution hyperspectral image, however, is not explicitly estimated. In this paper, we describe the application of a maximum a posteriori (MAP) estimation framework [14] with a stochastic mixing model (SMM) to the hyperspectral resolution enhancement problem. The stochastic mixing model provides a physical model to constrain the estimation problem in a similar manner as the linear mixing model in the approaches described in [10]–[12]. However, by using a statistical, as opposed to deterministic, framework, the approach does not rely on strict radiometric conformance and should be less sensitive to sensor degradations (noise, miscalibration, etc.) that reduce the radiometric correspondence between the two observations. The remainder of the paper is organized as follows. In Section II, the mathematical formulation of the MAP/SMM estimation approach is described, including comparison to the PCS and LSE methods (which will serve as performance benchmarks). In Section III, the characteristics of the test imagery and the quantitative metrics used to assess resolution enhancement performance are given. Experimental results are presented in Section IV, and conclusions are given in Section V.
of hyperspectral image pixels. This is a band-interleaved-bypixel spatially rastered representation of the inherently threedimensional (3-D) data. The objective of hyperspectral resolution enhancement is to estimate (or reconstruct) what is termed the high-resolution hyperspectral image, i.e., one possessing the spectral content of the low-resolution hyperspectral image but at the spatial resolution of the panchromatic image. This image is denoted as
II. MATHEMATICAL FORMULATION Suppose that an area is remotely imaged concurrently with a hyperspectral imaging system and a higher resolution panchromatic imaging system. The hyperspectral image observation will be termed “low-resolution” from the standpoint that its spatial resolution is lower than the panchromatic image. The panchromatic image is described as a one-dimensional (1-D) vector of intensity measurements (1) is the panchromatic image intensity at the spatial where is the number of location designated by the index , and panchromatic image pixels. This 1-D, or lexicographic, format is a rastered representation of the two-dimensional image. In a similar manner, the low-resolution hyperspectral image is also represented as a 1-D vector of spectroradiometric measurements
(3) where is a high-resolution pixel spectrum (with bands) at the spatial location designated by index , and is the number of image pixels (both corresponding to the panchromatic image). This is again a band-interleaved-by-pixel spatially rastered representation of the inherently 3-D data. A. Observation Model The relationship between the panchromatic image and the high-resolution hyperspectral image can be written in matrix form as (4) where the spectral response matrix is a sparse matrix whose columns are the spectral response functions for the panchromatic pixel locations, and is assumed to be a spatially independent zero-mean Gaussian random process with a standard deviation of . Similarly, the low-resolution hyperspectral image is related to the high-resolution hyperspectral image by (5) where is a sparse matrix whose rows are the spatial response functions for the low-resolution hyperspectral pixels and is assumed to be a spatially independent zero-mean Gaussian . Note random process with a spectral covariance matrix that both and are noninvertible, nonsquare matrices, so there is no unique solution for in either equation even in the absence of observation noise. B. MAP Formulation The goal of MAP estimation, as applied to the hyperspectral resolution problem, is to find an estimate for the high-resolution hyperspectral image that maximizes its conditional probability relative to the two observations, or (6) Using Bayes’ rule, the conditional probability density function can be expressed as (7)
(2) where is the low-resolution pixel spectrum (with bands) at the spatial location designated by index , and is the number
If the random noise processes of the two observations are not correlated (generally a good assumption given that the panchromatic and hyperspectral data come from different sensors), then
1926
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
the observed images are conditionally independent and can be written as the product of two conditional density functions (8)
Note that this function is quadratic in , so that a closed-form solution for the minimum can be readily obtained by determining the where the gradient of the cost function is zero, and its second derivative is positive. With some mathematical manipulation, it can be shown that this estimate corresponds to the solution to the linear system of equations
Ignoring the term in the denominator, which does not depend on , the MAP estimate is then given by (16) (9) The conditional density functions are easily derived from the observation models discussed in Section I. They take the form of Gaussian density functions by virtue of the assumed noise statistics
where (17) and (18)
(10) and
The Hessian matrix is a square, positive semidefinite matrix. Therefore, a solution to (16) will always represent a minimum. D. Stochastic Mixing Model
(11) The density function for the high-resolution hyperspectral image is unknown. Assume, however, that it is normally distributed. Then, it can be expressed as
(12) where the mean vector and covariance matrix can, in general, vary as a function of spatial location. Assuming spatial indepenwill be block diagonal. More specifically, the mean dence, vector and covariance matrix will take the forms (13) and
The solution of (16) requires some underlying model for the high-resolution hyperspectral image, or more specifically, the mean vector and covariance matrix as a function of spatial position. The stochastic mixing model is derived from the low-resolution hyperspectral image and used for this purpose. The formulation and implementation of the stochastic mixing model is described in this subsection. For a more detailed description, see the basic formulation in [16] and refinements in [17]. The stochastic mixing model is similar to the well-developed linear mixing model in that it attempts to decompose spectral data in terms of a linear combination of endmember spectra. However, in the case of the stochastic mixing model, the endare normally distributed random vectors, parammembers and covariance matrix eterized by their mean vector , instead of deterministic spectra. The endmember index ranges from one to , the number of endmembers assumed of the observed hyperspectral in the scene. Each spectrum image is then assumed to be a member of a mixture class defined by the random vector , where
(14) where
(19)
refers to a direct sum [15].
C. Optimization From the form of the probability density functions that compose (9), an objective function can be formed whose minimization will result in the optimal estimate. This objective function is given by the exponents of (10)–(12)
(15)
are the endmember abundances associated with the mixture class, and is the mixture class index. The endmember abundances are assumed to conform to the physical constraints of being between zero and one, and summing to unity, and are quantized at a specified level. This is done by a structured, recursive search algorithm that determines all the possible combinations of quantized abundances (i.e., discrete levels between zero and one) for the specified number of endmembers that satisfy the physical constraints. The result is a finite and discrete set of mixture classes.
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
1927
Due to the linear relationship in (19), the mean vector and covariance matrix for the th mixture class are functionally dependent on the corresponding endmember statistics by (20) and (21) and the class-conditional probability density function is described by the Gaussian mixture distribution Fig. 1. Graphical depiction of conventional LSE method.
(22) where
(23) The mixture model parameters that must be estimated from the measured data are the prior probabilities for all the mixture classes and the mean vectors and covariance matrices of the endmember classes. The discrete set of feasible classes are prede. fined by the set of abundances The stochastic expectation maximization [18] approach is used to self-consistently estimate the model parameters , , and , and the mixture class assignments [i.e., abundance maps ] for the observed low-resolution hyperspectral scene data. The specific implementation used in this paper incorporates the modifications detailed in [17]. It proceeds in the following manner. First, the endmember statistics are initialized. The mean vectors are initialized to a set of well-separated image spectra selected from the image by simplex maximization [19], the covariance matrices are initialized to the global covariance matrix, and the prior probabilities are uniformly distributed. Then, a series of iterations (nominally 20) are performed, each consisting of the following three steps. Step 1) The posterior probability is estimated on the th iteration for all combinations of image spectra and mixture classes according to (24)
Step 2) The image spectra are classified into the mixture classes in (19) by a Monte Carlo process based on the posterior probability values. Step 3) The endmember mean and covariance estimates are updated based on the sample statistics of spectra classified into the respective pure classes, and the prior probability estimates are updated based on the relative populations of the mixture classes.
After the algorithm converges, the low-resolution image spectra are ultimately assigned to the mixture classes , , that maximize the conditional posterior probability in (24), and low-resolution abundance maps are formed. That is, there are abundance maps giving the spatial abundance distributions for each endmember. The estimated abundance maps are then bilinearly , interpolated to give the high-resolution abundance maps , . According to (20) and (21), the elements of the spatially varying statistics in (13) and (14) are then estimated by (25) and (26) This provides the remaining terms in (17) and (18) to solve the system of equations in (16) for the high-resolution hyperspectral image estimate. E. Comparison to LSE In this paper, an LSE solution to the resolution enhancement problem is used as a performance benchmark. The LSE estimate takes the form (27) is the average of those panchromatic intensity values where that span the portion of the scene corresponding to . The vector can be estimated as a cross covariance [6] or using a least squares linear regression technique [1] based on and the available lower resolution images (e.g., the set of samples). The cross-covariance approach is corresponding used in this paper. A graphical comparison of the LSE and MAP/SMM methods is given in Figs. 1 and 2 for a simple case consisting of a two and corresponding to a panchromatic intensity values single, two-band, low-resolution spectrum . The LSE method constrains the estimated high-resolution spectra and to the
1928
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
III. DATA AND METRICS A. Test Data
Fig. 2.
Graphical depiction of proposed MAP/SMM method.
indicated line given by (27). The estimates are symmetrically located about (i.e., their average spectrum equals ) with a separation proportional to the difference in panchromatic intensity values. The MAP/SMM method statistically balances consistency with the observation models for and along with the mixture class probability density function from the interpolated abundance maps and SMM endmember statistics. This is depicted in Fig. 2 for the same simple case and uniform spatial and spectral response functions. Consistency with the model for forces the estimated high-resolution spectra toward the two parallel lines in the figure while consistency with occurs when the estimates are symmetrically located about . The SMM forces the estimates into their respective mixture classes. Variation can result toward the corresponding endmember classes and occurs particularly when neighboring pixels are purer (due to the linear interpolation of the abundance maps). F. Computational Issues The Hessian matrix is an matrix, so (14) represents a very large system of equations. It is also not necessarily a well-conditioned system of equations because of the high degree of correlation that often exists in hyperspectral data. That is, the rank of the covariance matrices associated with hyperspectral imagery is often less than the number of spectral bands. This issue is addressed by processing in the principal components domain rather than the spectral domain and limiting the data to a leading order subspace corresponding to the highest eigenvalues. Most hyperspectral imagery can be well represented with 10–20 principal components. In this paper, we also make the assumption that the point spread function corresponding to neighboring low-resolution pixels is nonoverlapping. Optically, this corresponds to the case where the resolution of the low-resolution hyperspectral image is strictly defined by system of equations the detector area. It allows the of smaller ( to be decomposed into a number ) systems equations, which significantly reduces the required computation. Computationally efficient methods to solve the problem without employing this assumption, however, are currently being studied.
The resolution enhancement results presented in this paper are based on three hyperspectral data sources. The first is a synthetically generated image, which provides a well-controlled test case where the imagery is known to conform to the basic model assumptions. The image was generated by the Digital Imaging and Remote Sensing Synthetic Image Generator (DIRSIG) model [20]. DIRSIG is a rigorous end-to-end spectroradiometric scene model that incorporates full 3-D scene construction; radiation propagation, absorption, and scattering; and extensive target-background interactions. The particular scene used is an emulation of data from the Hyperspectral Digital Imagery Collection Experiment (HYDICE) airborne sensor [21] viewing a forested area containing a clearing with a dirt road. The spatial resolution is approximately 0.15 m, about four to six times better than that collected by the HYDICE sensor during low-altitude operations (1.5–2-km altitude with 0.5-mrad angular resolution). The spectral range (0.38–2.6 m) and resolution (ranging from 4–20 nm) are based on the HYDICE sensor characteristics. The second image was collected by the Airborne Visible-Infrared Imaging Spectrometer (AVIRIS) sensor [22]. AVIRIS is a scanning dispersive hyperspectral imaging sensor that flies on the NASA ER-2 aircraft nominally 20 km above ground level and covers the visible/near-infrared (VNIR) and shortwave-infrared (SWIR) spectral regions. The spatial resolution is about 20 m, the spectral range is from 0.37–2.5 m, and the spectral resolution is on the order of 10 nm. The image used in this study was collected over Yorktown, VA. Because of the degree of complexity and low resolution of the image, it poses a challenge for resolution enhancement. The third image was collected by the NVIS airborne sensor [23]. NVIS is a pushbroom dispersive imaging sensor that flies on the Twin Otter aircraft nominally 1 km above ground level and also covers the VNIR/SWIR spectral range. The spatial resolution is about 0.6 m, the spectral range is from 0.45–2.5 m, and the spectral resolution is on the order of 5–10 nm. The image used in this study was collected over Fort A.P. Hill, VA, on May 4, 2000. The analysis presented in this paper is based only on the VNIR hyperspectral data covering the 0.45–0.9 m spectral range. In the DIRSIG and AVIRIS cases, the original imagery was spatially degraded by a factor of four to produce the low-resolution hyperspectral observation, and spectrally integrated over the entire spectral range to produce the high-resolution panchromatic observation. This produces imagery representative of hypothetical sensors with poorer spatial resolution and higher SNR (due to the spatial averaging) relative to the original imagery. This means that results obtained using this approach are not representative of data from the AVIRIS sensor were it to be used directly for the low-resolution hyperspectral data. Thus, the quantitative conclusions drawn are limited to comparing the MAP/SMM and LSE methods, and not the achievable enhancement for realistic AVIRIS data. The NVIS sensor also contains a bore-sighted VNIR panchromatic sensor of higher spatial resolution (0.5 ft), making it a good test set for the realistic employment of the resolution
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
Fig. 3. Panchromatic image for DIRSIG test case.
Fig. 4. Panchromatic image for AVIRIS test case.
1929
spatial and spectral characteristics, the metrics used to quantify the performance are designed to capture the difference in these fine details. Quantitative analysis of the NVIS test case is not possible because the lack of a truth image. The rms difference between the high-resolution hyperspectral image estimate and the truth image is the most basic error metric used. It is computed on a band-by-band basis. The spectral angle difference between the estimate and the truth for each spatial position is also computed. This metric isolates the errors in the spectral shape of the estimate from (correlated) intensity errors. In this paper, histograms of the spectral angle difference are used to compare methods. Neither of these metrics directly measures the ability of the hyperspectral resolution enhancement methods to reconstruct spatial content in the imagery beyond the resolution of the lowresolution hyperspectral sensor, the primary objective of the algorithms. To quantify this ability, a metric was developed that measures the correlation of the subpixel image structure to the truth. This is performed in the following manner. First, the estimate and the truth image are transformed to the principal component domain using the eigenvector matrix associated with the low-resolution spectral covariance matrix. Next, the mean intensity for each set of spatial pixels associated with a low-resolution pixel (i.e., super-pixel) is computed and subtracted. This is performed independently for each principal component band. Finally, the subpixel correlation metric is computed on a band-byband basis as the ratio of the cross correlation between the highpass filtered estimate with the truth image relative to the autocorrelation of the truth image
(28)
Fig. 5. Panchromatic image for NVIS test case.
enhancement method. Spatial registration of the panchromatic to hyperspectral imagery was performed by a three-step process consisting of local phase correlation measurement of relative image misregistration, spline interpolation of the measurements, and resampling with a finite-impulse-response interpolation kernel [24], [25]. The panchromatic imagery for the three datasets are illustrated in Figs. 3–5. B. Performance Metrics Four metrics were developed to compare the estimate produced by the resolution enhancement methods to the high-resolution truth image. Because the difference between the resolution enhancement methods lies in subtle, but important, fine
is the filtered estimate for the th spatial position and where the th principal component, and is the filtered truth for the same position and principal component. In the case of the DIRSIG forest scene, a spectral matched filter was constructed to detect the dirt road, a fine-resolution image feature that serves as a fourth test of the resolution enhancement methods. The filter vector was computed as the mean spectrum over the first 50 pixels in column 164, aligned directly on a soil track in the dirt road. Because this image is dominated by vegetation, which has fairly homogeneous spectral content, the spectral matched filter works very well at detecting the soil feature in the high-resolution image. The separation of the two soil tracks, however, is not resolved in the low-resolution result. Receiver operating characteristic (ROC) curves were generated in the following manner. First, a truth mask was made by thresholding the spectral matched filter output for the true high-resolution image at a level that visibly captured the soil track feature. A matched filter output was generated for each high-resolution image estimate and thresholded at a range of levels. The outputs above threshold at each level were then compared to the truth mask and scored in terms of the fraction detected and false-alarm rate to form the ROC curve.
1930
Fig. 6.
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
RMS difference for DIRSIG test case.
Fig. 7. Spectral angle error for DIRSIG test case.
IV. RESULTS A. DIRSIG Imagery Resolution-enhanced hyperspectral imagery using the PCS, LSE, and MAP/SMM methods was produced and quantitatively analyzed using the aforementioned metrics. These results are given in Figs. 6–9, where the SMM consisted of four endmem) and 74 mixture classes (i.e., ). There bers (i.e., is a modest improvement in the rms error relative to LSE with the MAP/SMM method, and a substantial improvement relative to PCS. Note that wavelet-based merging of the panchromatic image into the first principal component [26] is likely to produce a result similar to LSE. Similar trends occur for the spectral angle error results, where the MAP/SMM method provides a significant reduction in the number of pixels with spectral angle error in the 10 to 30 range relative to LSE. The subpixel cross-correlation results support the hypothesis that incorporation of the SMM will support the accurate reconstruction of subpixel spectral content, even down to the lower principal components. Enhancement for the LSE method is limited to the first
Fig. 8. Subpixel cross correlation for DIRSIG test case.
Fig. 9.
ROC curve for soil track detectability for DIRSIG test case.
principal component. Because detection of the soil track feature relies on high-resolution spectral information, the ROC curve associated with this detection problem is significantly improved by the MAP/SMM resolution enhancement. This is probably the best indicator of the efficacy of the MAP/SMM approach for this dataset. Qualitatively, the MAP/SMM method was found to be able to reconstruct subpixel spatial structure in the lower order principal components, while the LSE method is only able to enhance the first principal component. There was no visible difference between the first principal component imagery for both cases, but there were substantial differences in the lower order components. For example, Fig. 10 compares the second principal components for both cases as compared to the low-resolution observation and high-resolution truth. B. AVIRIS Imagery The results based on the AVIRIS imagery were fairly consistent with the DIRSIG results. Again, the low-resolution SMM
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
(a)
1931
(b)
Fig. 12. Subpixel cross correlation for AVIRIS test case.
(c)
(d)
Fig. 10. Second principal component images for DIRSIG test case. (a) Low resolution. (b) LSE estimate. (c) MAP/SMM estimate. (d) High-resolution truth.
Fig. 11.
Spectral angle error for AVIRIS test case.
(a)
(b)
(c)
(d)
Fig. 13. Third principal component images for AVIRIS test case. (a) Low resolution. (b) LSE estimate. (c) MAP/SMM estimate. (d) High-resolution truth.
consisted of four endmembers and 74 mixture classes. A subset of the quantitative and qualitative results is given in Figs. 11–13. Similarly to the DIRSIG case, the MAP/SMM method was able to reconstruct fine-resolution structure in the lower order principal components that LSE was not. This is clearly evident in the magnified, third principal component images in Fig. 13. Relative to the DIRSIG results, the quantitative improvement in spectral angle error and subpixel cross correlation compared to LSE is more modest for the AVIRIS test case.
only the first principal component contained significantly enhanced spatial resolution relative to the low-resolution image (not shown). The lower principal components were visibly identical to the low-resolution image. There is, however, higher resolution content in all principal components images for the MAP/SMM result. In this case, there was no high-resolution image truth to provide a quantitative comparison. Instead, the purpose was to demonstrate the efficacy of the algorithm in the real situation, i.e., coincident imagery from separate sensors.
C. NVIS Imagery
D. Sensitivity Analysis
Again, four endmembers and 74 mixture classes were used to derive the SMM for the NVIS case. The resulting principal component images are compared in Fig. 14. For the LSE result,
A sensitivity analysis was performed for both the DIRSIG and AVIRIS test cases to assess the sensitivity of MAP/SMM resolution enhancement performance to panchromatic image noise
1932
IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 42, NO. 9, SEPTEMBER 2004
(a)
Fig. 15. Soil track detectability versus panchromatic image SNR for DIRSIG test case.
(b)
(c)
Fig. 16. Soil track detectability versus hyperspectral image resolution factor for DIRSIG test case.
(d) Fig. 14. Principal component images for NVIS test case. (Left) LSE and (right) MAP/SMM. (a) First principal component. (b) Second principal component. (c) Third principal component. (d) Fourth principal component.
and hyperspectral image spatial resolution. The results are presented here in terms of the DIRSIG test case and the soil track detectability metric (the most sensitive metric), but the trends for the AVIRIS imagery and the other metrics were consistent. Fig. 15 shows the detectability results for a range of panchromatic image SNRs, defined as the ratio of the image-to-noise standard deviation. There is no reduction in performance until
the SNR drops into the 10–100 range, but the detectability drops quickly below an SNR of 10. Fig. 16 illustrates the reduction in performance as the spatial resolution of the hyperspectral observation image is varied while the panchromatic image is not changed. The legend indicates the ratio of panchromatic to hyperspectral spatial resolutions. The baseline case is a resolution factor of 4. A steady decrease in detectability with hyperspectral image spatial resolution is observed. V. CONCLUSION The results presented in this paper have provided an experimental proof-of-concept that the incorporation of a stochastic mixing model derived from the low-resolution hyperspectral observation in a maximum a posteriori estimator can provide better hyperspectral resolution enhancement results than conventional methods such as principal component substitution and
EISMANN AND HARDIE: APPLICATION OF THE STOCHASTIC MIXING MODEL
1933
least squares estimation. For all of the datasets investigated, the MAP/SMM method was able to produce more accurate subpixel structure in the lower order principal component images than the conventional methods. Additionally, the method was found to be robust to variations (e.g., noise and spatial resolution) in the observation imagery. We anticipate that such resolution enhancement will have positive impacts on subsequent hyperspectral processing and analysis, especially where finer spatial resolution is desired. This, however, requires further investigation. Other areas where further work is warranted include generalizing the implementation to the case of overlapping point spread functions, incorporating multispectral high-resolution imagery, and thoroughly analyzing the noise characteristics of the enhanced estimate relative to the low-resolution observation. We expect some inherent reduction in SNR with resolution enhancement, but have not yet characterized this effect.
[14] R. C. Hardie and M. T. Eismann, “MAP estimation for cross-sensor spatial resolution enhancement,” IEEE Trans. Image Processing, vol. 13, pp. 1174–1184, Sept. 2004. [15] H. Eves, Elementary Matrix Theory. New York: Dover, 1966, p. 107. [16] A. D. Stocker and A. P. Schaum, “Application of stochastic mixing models to hyperspectral detection problems,” Proc. SPIE, vol. 3071, pp. 47–60, 1997. [17] M. T. Eismann and R. C. Hardie, “Initialization and convergence of the stochastic mixing model,” Proc. SPIE, vol. 5159, pp. 307–318, 2003. [18] P. Masson and W. Pieczynski, “SEM algorithm and unsupervised statistical segmentation of satellite images,” IEEE Trans. Geosci. Remote Sensing, vol. 31, pp. 618–633, May 1993. [19] M. E. Winter, “Fast autonomous endmember determination in hyperspectral data,” in Proc. 13th Int. Conf. Applied Geological Remote Sensing, vol. II, 1999, pp. 337–344. [20] J. R. Schott et al., “Incorporation of time-dependent thermodynamic model and radiation propagation model into infrared three-dimensional synthetic image generation,” Opt. Eng., vol. 37, no. 7, p. 1505, July 1992. [21] R. R. Basedow et al., “The HYDICE instrument design,” in Proc. ISSR, Nov. 1992. [22] R. O. Green et al., “Imaging spectroscopy and the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS),” Remote Sens. Environ., vol. 65, pp. 227–248, 1998. [23] C. Simi et al., “Night Vision Imaging Spectrometer (NVIS) performance parameters and their impact on various detection algorithms,” Proc. SPIE, vol. 4049, pp. 218–229, 2000. [24] A. Schaum and M. McHugh, “Analytic methods of image registration: Displacement estimation and resampling,” Nav. Res. Lab., West Bethesda, MD, NRL Rep. 9298, 1991. [25] P. Landcaster and K. Salkauskas, Curve and Surface Fitting: An Introduction. London, U.K.: Academic, 1986. [26] T. Ranchin and L. Wald, “Fusion of high spatial and spectral resolution images: the ARSIS concept and its implementation,” Photogram. Eng. Remote Sens., vol. 66, pp. 49–56, 2000.
ACKNOWLEDGMENT The authors would like to acknowledge the support of R. Bartell (Veridian Corporation), B. Kendall (Space Computer Corporation), E. Winter (Technical Research Associates), and the NASA Jet Propulsion Laboratory for the imagery used in this paper. REFERENCES [1] J. C. Price, “Combining panchromatic and multispectral imagery from dual resolution satellite instruments,” Remote Sens. Environ., vol. 21, pp. 119–128, 1987. [2] C. K. Munechika, J. S. Warnick, C. Salvaggio, and J. R. Schott, “Resolution enhancement of multispectral image data to improve classification accuracy,” Photogram. Eng. Remote Sens., vol. 59, pp. 67–72, 1993. [3] W. J. Carper, T. M. Lillesand, and R. W. Kiefer, “The use of intensity-hue-saturation transformations for merging SPOT panchromatic and multispectral image data,” Photogram. Eng. Remote Sens., vol. 56, pp. 459–467, 1990. [4] V. K. Shettigara, “A generalized component substitution technique for spatial enhancement of multispectral images using a higher resolution data set,” Photogram. Eng. Remote Sens., vol. 58, pp. 561–567, 1992. [5] A. E. Iverson and J. R. Lersch, “Adaptive image sharpening using multiresolution representations,” Proc. SPIE, vol. 2231, pp. 72–83, 1994. [6] R. Nishii, S. Kusanobu, and S. Tanaka, “Enhancement of low resolution image based on high resolution bands,” IEEE Trans. Geosci. Remote Sensing, vol. 34, pp. 1151–1158, Sept. 1996. [7] D. P. Filberti, S. E. Marsh, and R. A. Schowengardt, “Synthesis of imagery from high spatial and spectral resolution from multiple image sources,” Opt. Eng., vol. 33, no. 8, Aug. 1994. [8] B. Zhukov, M. Berger, F. Lanzl, and H. Kaufman, “A new technique for merging multispectral and panchromatic images revealing sub-pixel spectral variation,” in Proc. IGARSS, Firenze, Italy, July 10–14, 1995, pp. 2154–2156. [9] B. Zhukov, D. Oertel, and F. Lanzl, “A multiresolution, multisensor technique for satellite remote sensing,” in Proc. IGARSS, Firenze, Italy, July 10–14, 1995, pp. 51–53. [10] H. H. Gross and J. R. Schott, “Application of spectral mixture analysis and image fusion techniques for image sharpening,” Remote Sens. Environ., vol. 63, pp. 85–94, 1998. [11] G. D. Robinson, H. N. Gross, and J. R. Schott, “Evaluation of two applications of spectral mixing models to image fusion,” Remote Sens. Environ., vol. 71, pp. 272–281, 2000. [12] H. N. Gross and J. R. Schott, “Application of spatial resolution enhancement and spectral mixture analysis to hyperspectral images,” Proc. SPIE, vol. 2821, pp. 30–41, 1996. [13] M. E. Winter and E. M. Winter, “Physics-based resolution enhancement of hyperspectral data,” Proc. SPIE, vol. 4725, pp. 580–587, 2002.
Michael T. Eismann (M’03) received the B.S. degree in physics from Thomas More College, Crestview Hills, KY, in 1985, the M.S. degree in electrical engineering from the Georgia Institute of Technology, Atlanta, in 1987, and the Ph.D. degree in electrooptics from the University of Dayton, Dayton, OH, in 2004. He is the currently a Technical Advisor with the Electro-Optical Targeting Branch, Sensors Directorate, Air Force Research Laboratory (AFRL), WrightPatterson Air Force Base, OH. He is responsible for overseeing the development of passive electrooptical and infrared sensor technology, and transition into operational airborne targeting and reconnaissance systems. His recent research emphasis has been the development of hyperspectral imaging sensors and image processing. Prior to joining AFRL in 1996, he was with the Environmental Research Institute of Michigan (ERIM). He joined ERIM in 1991, where he was involved in research concerning active and passive optical and infrared targeting and reconnaissance, optical information processing, and holographic optics.
Russell C. Hardie (S’91–M’92–SM’99) received the B.S. degree (magna cum laude) in engineering science from Loyola College, Baltimore, MD, in 1988, and the M.S. and Ph.D. degrees in electrical engineering from the University of Delaware, Newark, in 1990 and 1992, respectively. He is currently an Associate Professor with the Department of Electrical and Computer Engineering, University of Dayton, OH, and holds a joint appointment with the Electro-Optics Program. He was a Senior Scientist with the Earth Satellite Corporation prior to his appointment at the University of Dayton in 1993. His research interests include a wide variety of topics in the area of digital signal and image processing, including image restoration and medical image processing. Dr. Hardie was a corecipient of the Rudolf Kingslake Medal and Prize from SPIE in 1998, for work on multiframe image resolution enhancement algorithms. In 1999, he received the School of Engineering Award of Excellence in Teaching at the University of Dayton. He received the First Annual Professor of the Year Award in 2002 from the student chapter of the IEEE at the University of Dayton.
