Potential of Resolution-Enhanced Hyperspectral Data for Mineral Mapping Using Simulated EnMAP and Sentinel-2 Images

Abstract

Spaceborne hyperspectral images are useful for large scale mineral mapping. Acquired at a ground sampling distance (GSD) of 30 m, the Environmental Mapping and Analysis Program (EnMAP) will be capable of putting many issues related to environment monitoring and resource exploration in perspective with measurements in the spectral range between 420 and 2450 nm. However, a higher spatial resolution is preferable for many applications. This paper investigates the potential of fusion-based resolution enhancement of hyperspectral data for mineral mapping. A pair of EnMAP and Sentinel-2 images is generated from a HyMap scene over a mining area. The simulation is based on well-established sensor end-to-end simulation tools. The EnMAP image is fused with Sentinel-2 10-m-GSD bands using a matrix factorization method to obtain resolution-enhanced EnMAP data at a 10 m GSD. Quality assessments of the enhanced data are conducted using quantitative measures and continuum removal and both show that high spectral and spatial fidelity are maintained. Finally, the results of spectral unmixing are compared with those expected from high-resolution hyperspectral data at a 10 m GSD. The comparison demonstrates high resemblance and shows the great potential of the resolution enhancement method for EnMAP type data in mineral mapping.

1. Introduction

Hyperspectral images provide important information for the characterization of lithology. Using hundreds of contiguous spectral bands with a narrow bandwidth, a pixel spectrum is generated for the potential discrimination of target minerals [1,2]. Spectral unmixing is an important technique for discerning minerals at a subpixel scale using hyperspectral image data. Geological mapping and mineral exploration with hyperspectral images have been studied for decades [3,4,5]. In particular, the shortwave infrared (SWIR) spectral range of 2.0–2.5 µm is considered to have the most important diagnostic absorption features, which are related to the vibration processes in atomic bonds. Numerous distinctive minerals associated with primary lithologies and with various alteration styles can be identified; for instance, Al–OH at 2.2 µm, Mg–OH at 2.3 µm, and Ca–CO₃ at 2.32–2.35 µm [6].

Most geological applications using hyperspectral remote sensing, including mineral explorations, have relied on airborne acquisition to obtain high-spatial-resolution data. However, the airborne acquisition of hyperspectral images with high spatial resolution can become overly expensive when the coverage is very large. The spaceborne Environmental Mapping and Analysis Program (EnMAP) is dedicated to environmental applications of state-of-the-art hyperspectral technology (www.enmap.org) [7]. The launch of EnMAP is planned for mid-2018 and the expected mission lifetime is five years. A major objective of EnMAP is the monitoring and exploration of natural resources. Its fine spectral resolution of 6.5 nm in the visible and near infrared (VNIR) range (420–1000 nm) and 10 nm in the SWIR range (900–2450 nm) will make it well suited for geological mapping. It will have a swath width of 30 km and 242 spectral bands at 0.4–2.5 µm. This will be a great improvement on the 7 km swath width of the concurrent Hyperion [8]. EnMAP has a 27-day revisit cycle (off-nadir four days) and is able to acquire an along-track length of 5000 km per day based on the acquisition requests of users. Even equipped with the latest technology, a compromise must be made between scientific ambition and practicality regarding the spatial resolution. Therefore, EnMAP images will be acquired at a ground sampling distance (GSD) of 30 m. This spatial resolution will unfortunately limit the applicability of such data to a certain extent.

Owing to the high potential of spaceborne hyperspectral images in many applications, researchers are keen to look into possibilities for enhancing the spatial resolution. Some have studied superresolution image reconstruction using temporal and angular images [9,10]. The challenges in such approaches are the availability of quality temporal images and the complex manipulation of big data with high dimensionality. Changes in surface cover that occur between acquisitions can also have negative effects on the operation. Recent developments in image fusion have motivated new exploration of the spatial resolution enhancement of hyperspectral images [11,12]. Pansharpening, a traditional method applied to multispectral images, has received increasing attention for its use with hyperspectral images [13,14]. Its popularity can be attributed to the convenience of obtaining high-resolution multispectral images in the VNIR spectral range by a conventional method. The expected increasing demand for high-resolution hyperspectral images in the future has encouraged more studies, and recently newer resolution enhancement approaches based on hyperspectral and multispectral image fusion have shown highly promising results [15,16,17,18,19,20,21,22,23,24].

Sentinel-2 will provide high-resolution multispectral data with a short revisiting time (every 5 days with two satellites). The mission will acquire 13 bands in the VNIR and SWIR ranges with various spatial resolutions (10, 20 and 60 m), which contribute to a wide range of application including mineral mapping [25].

Table 1. Spectral bands for the Sentinel-2 sensor.

Band Number	Central Wavelength (nm)	Bandwidth (nm)	GSD (m)
1	443	20	60
2	490	65	10
3	560	35	10
4	665	30	10
5	705	15	20
6	740	15	20
7	783	20	20
8	842	115	10
8b	865	20	20
9	945	20	60
10	1380	30	60
11	1610	90	20
12	2190	180	20

lists the band descriptions. For experiments on resolution enhancement with hyperspectral EnMAP data, Sentinel-2 data is a very good option for obtaining the required high-resolution images. Even though changes between acquisitions still could occur, the 5-day revisiting time means there is a higher probability of avoiding such changes and it is easier to obtain high-resolution VNIR images at a similar time to hyperspectral acquisition.

A key limitation of previous studies on hyperspectral and multispectral image fusion is that fused images are validated mainly by quantitative measures of image quality but studies on its impact on applications are still lacking. An accurate validation of the fusion process on a pixel level can be performed only within a simulation study, because in such case the optimum result is known. Many of the conventional simulation studies consider mainly the image generation process with spectral response functions (SRFs), point spread functions (PSFs) and noise for simulation. However, a more realistic simulation should takes into account the entire image acquisition and processing chain because the characteristics of real data are influenced by instrumental and environmental parameters that naturally also impact the quality of the fused data products.

In this study, a simulated EnMAP-Sentinel image set forms the basis to investigate whether spatially enhanced hyperspectral data obtained from the fusion of hyperspectral and multispectral images can retain spectral fidelity and consequently contribute to mineral mapping. Sensor end-to-end simulation tools [26,27,28,29] that consider the entire image acquisition and processing chains are used for simulating realistic EnMAP and Sentinel-2 images. We demonstrate the fusion of simulated EnMAP and Sentinel-2 VNIR images using various methods (one matrix factorization method and two pansharpening methods) and validate the quality of the fused data products in regard to the preservation of the original spectral information content. Among the hyperspectral and multispectral image fusion methods, we have focused on a matrix factorization method that requires relative sensor configurations, such as SRFs and PSFs, but produces good results with only a moderate computational burden [14]. Our contribution is threefold: we present (1) a validation of hyperspectral image enhancement by continuum removal and spectral unmixing in terms of its applicability in mineral mapping; (2) image fusion using sensor end-to-end simulation data; (3) hyperspectral and multispectral image fusion without the overlap of SRFs.

2. Hyperspectral and Multispectral Image Fusion

2.1. Related Work

Hyperspectral and multispectral image fusion enables the production of higher-spatial-resolution hyperspectral imagery while preserving high spectral accuracy, making it advantageous over pansharpening. The fusion problem is formulated as an ill-posed inverse problem. Hyperspectral and multispectral image fusion is still an open issue, and the methods presented in the literature for tackling this issue are limited in number compared to those proposed for pansharpening.

A wavelet-based technique that inherited the pansharpening algorithm was first proposed for hyperspectral and multispectral image fusion [30,31]. However, its performance strongly depended on the spectral sampling method, making it difficult to enhance the spatial resolution of all hyperspectral band images.

A Bayesian method based on a maximum a posteriori (MAP) estimation was a major breakthrough in terms of the capability to enhance the spatial resolution of all hyperspectral band images using higher-spatial-resolution data acquired by a panchromatic or multispectral imager [15,16]. This method adopted a stochastic mixing model (SMM) to estimate the underlying spectral scene characteristics and integrated them into a MAP-based fusion framework. The optimization was performed in the principal component domain rather than the spectral domain to decrease the computational burden. The MAP/SMM method outperformed those based on conventional pansharpening algorithms and has been further extended using wavelet transforms to improve noise resistance [17].

Methods based on a matrix factorization or spectral unmixing perspective have been attracting attention in recent years owing to their straightforward interpretation of the fusion process and state-of-the-art performance [18,19,20]. Coupled nonnegative matrix factorization (CNMF) formulates the image fusion problem as a coupled spectral unmixing problem [18]. A high-spectral-resolution endmember matrix and high-spatial-resolution abundance matrix are obtained from hyperspectral and multispectral images, respectively, by coupled spectral unmixing. The fused image can be reconstructed by the multiplication of the two underlying matrices. Sparse regularization on the abundance fractions or local image processing was shown to be effective for mitigating the ill-posedness of the inverse problem and improving the quality of fusion products [19,20].

When the two underlying matrices of the fused image are generalized as the subspace transform matrix and the representation coefficient matrix, the fusion problem is formulated within the Bayesian fusion model, which includes the MAP/SMM method. Recent advances in this approach can be characterized as the use of regularization in the problem by defining an appropriate prior distribution for the observed scene, such as patch-based sparse representation and vector-total-variation-based regularization of the spatial distribution of subspace coefficients [21,22].

Pansharpening methods can be applied to hyperspectral and multispectral image fusion by regarding each multispectral band and spectrally corresponding hyperspectral bands as panchromatic and multispectral images, respectively, in the conventional pansharpening framework. A generalization of pansharpening that can be adapted to hyperspectral and multispectral image fusion has recently been studied [23]. One of the state-of-the-art pansharpening methods based on patch-based sparse representation was successfully extended to solve the fusion problem [24]. Major challenges for pansharpening-based methods are the discontinuity of spectra in the fused image and the fusion of spectrally non-overlapping images.

Many recent methods are summarized in a recently published review paper of hyperspectral pansharpening [14], which can be seen as a particular case of hyperspectral and multispectral image fusion. In this paper, we focus on CNMF owing to its straightforward interpretation, robust performance as a benchmark method as shown in the literature, and its applicability to the fusion of spectrally non-overlapping images, which will be tackled in this paper.

2.2. The CNMF Algorithm

Hyperspectral and multispectral image fusion aims to estimate a high-spatial-resolution hyperspectral image ($\mathbf{X}\in {ℝ}^{L\times P}$) by fusing a low-spatial-resolution hyperspectral image (${\mathbf{Y}}_h\in {ℝ}^{L\times P_h}$) and a high-spatial-resolution multispectral image (${\mathbf{Y}}_m\in {ℝ}^{L_m\times P}$). L and $L_m$ denote the numbers of spectral bands of the hyperspectral and multispectral sensors, respectively. $P_h$ and P denote the numbers of pixels in the hyperspectral and multispectral images, respectively. $L>L_m$ and $P_h<P$ are satisfied by considering the trade-off between the spectral and spatial resolutions of the two sensors.

CNMF aims to reconstruct high-spatial-resolution hyperspectral data by extracting the underlying spectral scene characteristics at a high spatial resolution via matrix factorization, which can be interpreted as spectral unmixing. Spectral unmixing of the hyperspectral image provides the spectral signatures of endmembers and a set of fractional abundances for each pixel at a large GSD. The multispectral image itself cannot achieve accurate spectral unmixing due to broader SRFs at limited wavelength ranges; however, it is helpful to determine the detailed spatial distribution of abundances at a smaller GSD once the spectral signatures and the fractional abundances are obtained for each pixel at a large GSD by spectral unmixing of the hyperspectral image. The multiplication of the spectral signatures and their high-resolution abundance maps leads to the reconstruction of high-spatial-resolution hyperspectral data.

We assume that the two observed images are acquired over the same area under the same atmospheric and illumination conditions and are geometrically coregistered. The low-spatial-resolution hyperspectral image and the multispectral image can be viewed as degraded versions of the target image in the spatial and spectral domains, respectively. Therefore, ${\mathbf{Y}}_h$ and ${\mathbf{Y}}_m$ are modeled as

$${\mathbf{Y}}_h=\mathbf{\text{XBS}}+{\mathbf{N}}_h$$

(1)

$${\mathbf{Y}}_m=\mathbf{\text{RX}}+{\mathbf{N}}_m$$

(2)

Here, $\mathbf{B}\in {ℝ}^{P\times P}$ is a matrix modeling the PSF of the hyperspectral sensor, $\mathbf{S}\in {ℝ}^{P\times P_h}$ is the spatial downsampling matrix and $\mathbf{R}\in {ℝ}^{L_m\times L}$ is the spectral transform matrix, with its columns representing the spectral responses of the multispectral sensor. ${\mathbf{N}}_h$ and ${\mathbf{N}}_m$ are the residuals. The relative sensor characteristics, such as B and R, can be estimated from the observed data sources when they are not given as prior knowledge [32].

The unmixing-based approaches to fusion introduce linear mixture models for hyperspectral images. Linear spectral mixture models are commonly used for spectral unmixing owing to their physical effectiveness and mathematical simplicity [33]. The spectrum at each pixel is assumed to be a linear combination of several endmember spectra. Therefore, X is formulated as

$$\mathbf{X}=\mathbf{\text{EA}}+\mathbf{N}$$

(3)

where $\mathbf{E}\in {ℝ}^{L\times M}$ is the spectral signature matrix, with the column vectors corresponding to the M endmember spectra, $\mathbf{A}\in {ℝ}^{M\times P}$ is the abundance matrix, with each column vector denoting the abundance fractions of all the endmembers at the pixel, and N is the residual. The endmember spectra and abundance fractions are nonnegative. By substituting Equation (3) into Equations (1) and (2), ${\mathbf{Y}}_h$ and ${\mathbf{Y}}_m$ can be approximated as

$${\mathbf{Y}}_h\approx \mathbf{E}{\mathbf{A}}_h$$

(4)

$${\mathbf{Y}}_m\approx {\mathbf{E}}_m\mathbf{A}$$

(5)

where we define the spatially degraded abundance matrix ${\mathbf{A}}_h\in {ℝ}^{M\times P_h}$ and the spectrally degraded endmember matrix ${\mathbf{E}}_m\in {ℝ}^{L_m\times M}$ as

$${\mathbf{A}}_h=\mathbf{\text{ABS}}$$

(6)

$${\mathbf{E}}_m=\mathbf{\text{RE}}$$

(7)

CNMF focuses on the estimation of the spectral signatures of endmembers (E) and the corresponding high-resolution abundance maps (A) by the coupled spectral unmixing of the low-spatial-resolution hyperspectral image (${\mathbf{Y}}_h$) and the high-spatial-resolution multispectral image (${\mathbf{Y}}_m$). CNMF alternately unmixes ${\mathbf{Y}}_h$ and ${\mathbf{Y}}_m$ in the framework of nonnegative matrix factorization (NMF) [34,35] to estimate E and A under the constraints of the relative sensor characteristics given by Equations (6) and (7). NMF decomposes a nonnegative data matrix into a product of two nonnegative matrices, and it has been shown to be effective for spectral unmixing satisfying physical constraints without assuming the presence of pure pixels. The squared Frobenius norm of a residual matrix in a linear spectral mixture model is commonly used for a cost function. NMF unavoidably converges to local minima depending on the initialization. CNMF alternately takes advantage of the two images, i.e., the spectral resolution of the hyperspectral image and the spatial resolution of the multispectral image, to initialize the other spectral unmixing so that the optimization can converge to a better local minimum.

CNMF starts with NMF unmixing of the low-spatial-resolution hyperspectral image. The matrix E can be initialized using endmember extraction methods, for example, vertex component analysis (VCA) [36]. In this work, since our interest is focused on the SWIR range of fused data, this endmember extraction process is carried out using only SWIR bands. E and ${\mathbf{A}}_h$ are then alternately optimized by Lee and Seung’s multiplicative update rules [35]. Next, A is estimated from the multispectral image. ${\mathbf{E}}_m$ is set by Equation (7) and A is initialized by the spatially up-sampled matrix of ${\mathbf{A}}_h$ obtained by bilinear interpolation. The sequential unmixing for the hyperspectral image is performed after initializing ${\mathbf{A}}_h$ by Equation (6). After that, the two images are alternately unmixed until convergence. Finally, the target image is obtained by the multiplication of E and A. More details of the implementation are given in [18].

3. Materials and Validation Methods

3.1. Study Area and Data Preparation

An airborne HyMap scene of Rodalquilar in the south-eastern corner of Spain was used to generate EnMAP and Sentinel-2 Level-2a data products. The airborne image covers the Sierra del Cabo de Gata (Cabo de Gata National Park), where a gold mining area is situated that has been used in several studies [6,37]. Volcanic rocks have undergone a transformation to minerals such as silica, alunite, kaolinite, montmorillonite and chlorite. Various alteration phases can be distinguished, and gold deposits have been found in the central part of the volcanic field [38]. This area is commonly used to investigate mineral exploration applications. The top panel of Figure 1 illustrates the schematic surface alteration map of four types of hypogene (stage 1) alteration, i.e., silicic, advanced argillic, intermediate argillic, and propylitic, and one type of supergene (stage 2) acid sulfate alteration characterized by stage 2 alunite. The silicic alteration zone includes vuggy residual silica and massive silicified rock within halos of advanced argillically altered rocks. The advanced argillic alteration zone is composed mainly of quartz, (stage 1) alunite, and kaolinite. The intermediate argillic alteration zone covers a wide area of rocks with variable degrees of argillization and consists of quartz, kaolinite, illite, illite-smectite, K-feldspar, hematite, goethite, pyrite, and minor alunite, pyrophyllite, and diaspore. Argillically altered rocks grade into a propylitic assemblage that consists of quartz, chlorite, K-feldspar, vermiculite, illite, smectite, hematite, goethite, primary plagioclase and minor calcite. The bottom panel of Figure 1 shows a false-color composite image of the study area. Bright red represents cultivated cropland, and shades of red show distribution of various vegetated land covers. The vegetation is xerophytic and the landscape is characterized by shrubland.

The HyMap scene was acquired in June 2003. Owing to technical issues, the electromagnetic spectra at SWIR-1 (1.40–1.80 μm) are not available. The data was orthorectified with a 4 m GSD using PARGE and atmospherically corrected using ATCOR. On the basis of this dataset, EnMAP and Sentinel-2 L2a products (orthorectified surface reflectance data) were simulated using the sensor end-to-end simulation tools EeteS [26,27,28] and S2eteS [29], which were developed at the German Research Center for Geosciences (GFZ) in Potsdam, Germany. These tools simulate the entire image data acquisition, calibration and processing chain from spatially and spectrally oversampled data to intermediate Level-1a raw data and to the final EnMAP and Sentinel-2 products, such as Level-1b, Level-1c and Level-2a data. Figure 2 shows the flowchart of data preparation. The data acquisition consists of a sequential processing chain represented by four independent modules—atmospheric, spatial, spectral and radiometric. These modules allow flexible customization of a wide range of simulation input parameters. For example, an atmospheric simulation was performed for an image using typical atmospheric parameters: an aerosol optical thickness of 0.2, columnar water vapor of 1.5–2.5 g·cm⁻² and a rural aerosol model. The spatial and spectral modules include resampling of an image in the spatial and spectral domains using the sensor specific PSFs and SRFs, respectively. The radiometric module transformed the at-sensor radiance to DN by simulating instrumental noise and calibration coefficients. The modules are coupled with a backward simulation branch consisting of calibration modules such as nonlinearity, dark current and absolute radiometric calibration and a series of preprocessing modules such as radiometric calibration, co-registration, orthorectification and atmospheric correction. In case of Sentinel-2, all bands of the Level-1c- and Level-2a-products were spatially resampled to a uniform 10 m pixel grid. A reference hyperspectral image with a 10 m GSD for all bands was also prepared by applying only spatial and spectral resampling to the HyMap at-surface reflectance data. The HyMap data were spectrally resampled to a 1 nm resolution before the SRFs were applied for simulating the EnMAP and reference images.

While HyMap data have certain limitations pertaining to spectral and spatial resolutions, which is normal to any simulation procedure, the accurate simulation of EnMAP and Sentinel-2 data is possible with the robust parameter settings in our tools, which has been proven in previous reports [26,27,28,29]. Before real data are available, these well-tested simulation tools provide the closest resemblance to the realistic scenarios.

3.2. Validation Methods

The EnMAP and Sentinel-2 VNIR images were fused using the CNMF method. Owing to the popularity of the pansharpening approach for hyperspectral enhancement, we have included results from two representative pansharpening algorithms: (1) the Gram-Schmidt adaptive (GSA) algorithm [39]; and (2) MTF-tailored generalized Laplacian pyramid fusion (MTF-GLP) [40]. The most correlated Sentinel-2 VNIR image was used to sharpen each EnMAP SWIR image. Results from a single-frame spatial enhancement technique involving bicubic interpolation are also provided.

A reference image was used to validate the result of image fusion. We calculated the four quantitative measures used in [14] to check whether high image quality and spectral fidelity were maintained: the cross correlation (CC), spectral angle mapper (SAM), root mean square error (RMSE) and erreur relative globale adimensionnelle de synthèse (ERGAS) [41]. Let $\widehat{X}$ and X denote a fused image and a reference image, respectively.

(1)

CC is a characterization of geometric distortion obtained for each band with an ideal value of 1. We used an average value of CCs for all bands, which is defined as

$$\text{CC}\left(\widehat{X},\mathbf{X}\right)=\frac{1}{L}{\displaystyle \sum _{i=1}^L}\frac{{\sum }_{j=1}^P\left({\widehat{X}}_{ij}-\frac{1}{P}{\sum }_{j\prime =1}^P{\widehat{X}}_{ij\prime }\right)\left(X_{ij}-\frac{1}{P}{\sum }_{j\prime =1}^P{X}_{ij\prime }\right)}{\sqrt{{\sum }_{j=1}^P{\left({\widehat{X}}_{ij}-\frac{1}{P}{\sum }_{j\prime =1}^P{\widehat{X}}_{ij\prime }\right)}^2{\sum }_{j=1}^P{\left(X_{ij}-\frac{1}{P}{\sum }_{j\prime =1}^P{X}_{ij\prime }\right)}^2}}$$

(8)

(2)

SAM is a measure for the shape preservation of a spectrum calculated at each pixel with a unit degree and 0 as the ideal value. An average value of a whole image is defined as

$$\text{SAM}\left(\widehat{X},\mathbf{X}\right)=\frac{1}{P}{\displaystyle \sum _{i=1}^P}\text{arccos}\left(\frac{{\widehat{x}}_i^T{\mathbf{x}}_i}{\Vert {\widehat{x}}_i\Vert \Vert {\mathbf{x}}_i\Vert }\right)$$

(9)

where ${\mathbf{x}}_i$ denotes the $i^{th}$ column vector of X and ‖·‖ is the $l_2$ norm.

(3)

RMSE is calculated at each pixel as the difference of spectra between the fused image and the reference image. We used an average value of RMSEs for all pixels, which is defined as

$$\text{RMSE}\left(\widehat{X},\mathbf{X}\right)=\frac{1}{P}{\displaystyle {\sum }_{j=1}^P}\sqrt{\frac{1}{L}{\displaystyle {\sum }_{i=1}^L}{\left({\widehat{X}}_{ij}-X_{ij}\right)}^2}$$

(10)

For continuum-removed data, we consider a modified reference data ($X^{\prime }$) with each spectrum best matching its corresponding spectrum of the fused data to perform spectral feature fitting [42]. The modified continuum removed spectrum with scaling and shifting parameters at the $i^{th}$ pixel is obtained by the least squares algorithm as ${\mathbf{x}}_i^{\prime }=\left[{\mathbf{x}}_i 1_L\right]{({\left[{\mathbf{x}}_i 1_L\right]}^T\left[{\mathbf{x}}_i 1_L\right])}^{-1}{\left[{\mathbf{x}}_i 1_L\right]}^T{\mathbf{x}}_i$, where $1_L\in {ℝ}^{L\times 1}$ denotes a vector with all L entries one.

(4)

ERGAS provides a global statistical measure of the quality of fused data with the best value at 0, which is defined as

$$\text{ERGAS}\left(\widehat{X},\mathbf{X}\right)=100\sqrt{\frac{P}{P_h}}\sqrt{\frac{1}{L}{\sum }_{i=1}^L{\left(\sqrt{\frac{1}{P}{\sum }_{j=1}^P{\left({\widehat{X}}_{ij}-X_{ij}\right)}^2}/\frac{1}{P}\sum _{j=1}^P{X}_{ij}\right)}^2}$$

(11)

The quality of the fused data products was validated with reflectance and continuum removal data using these measures of image quality and visual interpretations. Continuum removal is a technique that highlights useful information such as the shape and location of absorption peaks from background noise. As mineral materials often have subtle absorption features, continuum removal has been successfully applied to reflectance spectra [43,44].

We later show that the unmixing performance of the enhanced image is comparable to that of a “real” high-resolution hyperspectral image. Since the ground-truth data at a 10 m pixel level were not available, it was not possible to evaluate the precision of the unmixing. Instead, and probably more important for our purpose, the unmixing results from the enhanced 10 m hyperspectral image were visually compared with those from the 10 m hyperspectral reference image. Pixel spectra are also presented to show the correctness of the image fusion method.

4. Results and Discussion

We fused Sentinel-2 bands 2, 3, 4 and 8 and EnMAP bands 1–127 and 189–242. In this way, a total of 181 high-resolution bands were simulated with 127 bands at 423–1335 nm, and 54 bands at 1998–2439 nm. Additionally, a 10 m hyperspectral reference image including the same EnMAP spectral bands was created. For mineral mapping, our analyses focused on the 52 SWIR bands between 2016–2439 nm. The quality indices are listed in

Table 2. Quality indices.

Data	Method	CC	SAM	RMSE	ERGAS
VNIR and SWIR	Cubic	0.91749	2.8365	0.01909	3.3857
	GSA	0.98629	2.713	0.01372	2.0112
	MTF-GLP	0.98571	2.692	0.01355	2.0151
	CNMF	0.988	2.6994	0.01349	1.9793
SWIR	Cubic	0.90549	1.6368	0.01694	2.9814
	GSA	0.97374	1.629	0.01149	1.9339
	MTF-GLP	0.97346	1.6381	0.01132	1.8958
	CNMF	0.97329	1.6193	0.01165	1.9336
Continuum removed SWIR	Cubic	0.76578	0.65553	0.01015	0.45248
	GSA	0.7745	0.64832	0.01011	0.4504
	MTF-GLP	0.7659	0.66141	0.01019	0.46186
	CNMF	0.83494	0.60761	0.00915	0.40899

in three categories: (1) average over all the fused bands (both VNIR and SWIR); (2) average over fused bands in only the SWIR range; and (3) average over fused bands in only the SWIR range with continuum removal. The two pansharpening methods and the CNMF method achieved very similar scores in categories 1 and 2, but in category 3 with continuum removal the latter significantly outperformed the pansharpening methods.

Figure 3 shows a false-color composite comparison of the enhanced images and the reference image, where typical absorption feature bands at 2259, 2201 and 2014 nm are represented by red, green and blue channels, respectively. Two close-up images are also shown. A visual comparison suggests high-quality reconstruction using both the pansharpening methods and the CNMF method, and their results are clearly superior to those of single-frame nearest-neighbor interpolation and bicubic interpolation. The presence of target minerals (e.g., alunite, kaolinite, and smectite), indicated by shades of pink, is enhanced with the detailed spatial structure recovered through the high-resolution Sentinel-2 VNIR bands. There are, however, no major differences among the different methods according to the visual comparison. Absolute reflectance (or radiance) data have been commonly used for visual comparison in the literature of the resolution enhancement of hyperspectral images. However, absolute reflectance values, which appear as brightness in Figure 3, are heavily influenced by topography and solar elevation angles rather than actual mineralogy.

Figure 4 shows the visuals after continuum removal to demonstrate spectral absorption features at 2201, 2159 and 2115 nm. Black, blue and light blue are corresponding to alunite, kaolinite and smectite, respectively. The low abundances of kaolinite are indicated by pale purple. These correspondences were drawn from the analysis of spectral unmixing, which will be shown later. It is clear that CNMF generated results with finer detail. This is consistent with and has been verified by the quality indices in Table 2. The color composite images of the GSA and MTF-GLP data are surprisingly as blurred as that of the bicubic interpolation image. This is because of the high spectral correlation of the injected spectral vectors in the GSA and MTF-GLP data. More precisely, the injected detail vectors are parallel to each other in the centered data in the case of GSA, and the injected detail vector at each pixel is parallel to the corresponding resampled hyperspectral vector in the case of MTF-GLP. Spatial details are represented by the magnitudes of the injected vectors, which are normalized by continuum removal, resulting in their disappearance. The fact that the two pansharpening methods have no effect on resolution after continuum removal (Figure 4), which is hidden in reflectance (Figure 3), implies that continuum removal is effective to validate the fused hyperspectral data products. To our knowledge, this is the first study to demonstrate the efficacy of continuum removal for assessing the added value of hyperspectral image enhancement and to reveal the advantage of matrix-factorization-based image fusion as well as limitations of pansharpening.

To demonstrate spectral differences between the reference and enhanced images, we selected some exemplary spectra. We used the information of the 21 points reported in [37] to identify regions with different mineralogy. Since the exact coordinates were not available, the positions of these 21 points were estimated as accurately as possible as shown in the enlarged image in the bottom panel of Figure 1. Figure 5 shows the reflectance spectra generated from the different methods. The endmember spectra from the reference image are smoother because the reference image was generated by spectral resampling of the HyMap at-surface reflectance image without atmospheric simulation and correction (Figure 2), whereas those of the enhanced image products include errors resulting from a non-optimal atmospheric correction of the EnMAP image that is performed in the EnMAP end-to-end simulation. The spectra generated by bicubic interpolation show greater deviation at several locations. Otherwise, very good matching is achieved by both the pansharpening methods and the CNMF method with the useful information preserved. In Figure 6, the spectral signatures after continuum removal are shown. Better matching with CNMF than with the other methods is apparent. This is consistent with the obtained quality indices (Table 2) and the visual interpretation of Figure 4. At 10 m, CNMF appears to have been very effective for image fusion, most importantly, the reconstruction captures the absorption features crucial for mineral identification.

To further illustrate the potential of the enhanced image, the results of a spectral mixture analysis were used. Assuming that the spectral signal at each pixel is a weighted combination of endmembers, spectral unmixing methods have been proposed to resolve the fraction or abundance of each individual endmember [45]. While a normal realistic scene contains a large number of components, each having a distinctive spectral signature, an image pixel normally has only a few of these components mixed together. Although each pixel is a mixture of different numbers of components, simple spectral unmixing methods use a fixed number of endmembers, which can be problematic. Therefore, multiple-endmember unmixing [46] has been proposed and applied to the Rodalquilar scene in [37].

The first step of our spectral mixture analysis was endmember extraction. We adopted VCA to extract endmember candidates from the reference data and a final set of endmembers were manually selected by visual inspection, which resulted in six classes: alunite, kaolinite, smectite, carbonate, nonphotosynthetic vegetation (NPV) and green vegetation (Figure 7). Since the location of the endmembers in [37] was unknown, the endmember extraction was carefully processed by referring the spectral signatures of endmembers and their abundance maps in [37]. We extracted 25 endmember candidates from the reference data using VCA and manually selected five endmembers (alunite, kaolinite, smectite, NPV and green vegetation) considering their spectral similarity to those presented in [37], their locations and the corresponding abundance fractions shown in [37]. Only the spectral signature of carbonate was extracted fully by manually checking its locations with a high abundance fraction given in [37]. The six extracted spectra were used to unmix the reference image. Spatially corresponding spectra in the enhanced data were used as endmembers for spectral unmixing of the enhanced data.

Next, multiple endmember spectral mixture analysis (MESMA) [46] was applied to each reflectance image using respective endmember spectra to estimate the abundance fractions. We used photometric shade, with zero reflectance in all channels, as the shade endmember. Forty-one models with two to four endmembers were formed: six models with two endmembers (signature + shade), fifteen models with three endmembers (two signatures + shade) and twenty models with four endmembers (three signatures + shade). The constraint of the abundance fractions was set as −0.01 ≤ f ≤ 1.01. The maximum allowable RMSE was set to 0.01 in reflectance. Color composite images of abundance fractions for alunite (red), kaolinite (green) and smectite (blue) obtained using the reference, CNMF, GSA, bicubic interpolation, and nearest-neighbor interpolation images are shown in Figure 8. Although our results were generated using a 10 m image, they show good indicative resemblance to the results obtained from the 4 m hyperspectral image in [37]. Our primary goal is to show that the abundance fractions obtained using the fused data are comparable to those obtained from the reference image. The abundance fractions for the target minerals generated from the CNMF image and the reference image are highly comparable, although some spatial details are missing in the results obtained from the CNMF image. On the other hand, the results obtained using the GSA and bicubic interpolation images are more blurred, including different distributions of alunite due to spectral distortion. These results show for the first time that small spectral distortions included in resolution-enhanced hyperspectral data may have a major impact on spectral unmixing. This visual validation of the results of spectral unmixing has confirmed the efficacy of the CNMF method in addition to validating the numerical quality of enhanced images. The good resemblance of the results of spectral unmixing with that reported in [37] indicates the potential of CNMF-fused data in mineral mapping applications.

5. Conclusions

In this work, we investigated the potential of resolution-enhanced hyperspectral data for mineral mapping using a simulated pair of EnMAP and Sentinel-2 images. The EnMAP and Sentinel-2 products were generated from a HyMap scene covering a mining area using sensor end-to-end simulation tools, which take into account various instrumental and environmental configurations. The EnMAP and four Sentinel-2 VNIR bands were fused to produce a 10 m GSD hyperspectral image using the matrix factorization method that reconstructs the fused image via extraction of common underlying spectral signatures and their corresponding abundance fractions from the two images. The quality of resolution-enhanced hyperspectral data was validated by comparison with state-of-the-art pansharpening methods. Numerical quality measures and the visual interpretations using continuum-removed data demonstrated the high spectral and spatial fidelity with the presented method. The results of spectral unmixing showed highly comparable distributions with those generated from the reference image at a finer spatial resolution. These results suggested the effectiveness and potential of the presented hyperspectral and multispectral image fusion method. The contribution is a synergetic use of future EnMAP and Sentinel-2 data resulted in a spatial definition in geological exploration finer than those of the current and near future satellite missions such as EnMAP, the Japanese Hyperspectral Imager Suite (HISUI) [47], the Italian PRISMA (Hyperspectral Precursor of the Application Mission) [48], the US Hyperspectral Infrared Imager (HyspIRI) [49] and the French HYPXIM [50]. Naturally, the image fusion technique presented in this work can be extended to conventional applications such as detailed land cover classification and target detection; and with the enhanced data at 10 m GSD, there are potentially new applications that have not been anticipated for spaceborne hyperspectral remote sensing before.

In this study, we only used the 10 m Sentinel-2 VNIR bands to minimize the complexity of scaling, resulting in a challenging image fusion task without the overlap of SRFs in the SWIR range. In future work, the incorporation of the six additional 20 m SWIR bands will be further investigated to mitigate the ill-posedness of the inverse problem and improve the quality of the resolution-enhanced products. In addition, further experiments will be conducted based on better simulated sensor images from APEX and HySpex data with higher spectral and spatial resolutions. Clearly, further research will be required to validate the image fusion method with multi-sensor real data. The Multi-User System for Earth Sensing (MUSES) platform on the international space station is expected to provide a pair of hyperspectral and multispectral images acquired over the same area under the same environmental conditions, which satisfies the assumption of image fusion [51].