HYPERSPECTRAL AND MULTISPECTRAL DATA FUSION BASED ON NONLINEAR UNMIXING Naoto Yokoya,1 Jocelyn Chanussot,2 and Akira Iwasaki3 1 Department of Aeronautics and Astronautics, The University of Tokyo, Japan GIPSA-lab, Signal and Image Dept., Grenoble Institute of Technology, Grenoble, France 3 Research Center for Advanced Science and Technology, The University of Tokyo, Japan 2 ABSTRACT Data fusion of low spatial-resolution hyperspectral (HS) and high spatial-resolution multispectral (MS) images based on a linear mixing model (LMM) enables the production of high spatial-resolution HS data with small spectral distortion. This paper extends the LMM based HS-MS data fusion to nonlinear mixing model using a bilinear mixing model (BMM), which considers second scattering of photons between two distinct materials. A generalized bilinear model (GBM) is able to deal with the underlying assumptions in the BMM. The GBM is applied to HS-MS data fusion to produce high-quality fused data regarding multiple scattering effect. Semi-nonnegative matrix factorization (Semi-NMF), which can be easily incorporated with the existing LMM based fusion method, is introduced as a new optimization method for the GBM unmixing. Comparing with the LMM based HSMS data fusion, the proposed method showed better results on synthetic datasets. Index Terms— Data fusion, nonlinear unmixing, bilinear mixing model, semi-nonnegative matrix factorization 1. INTRODUCTION Hyperspectral (HS) imaging sensors generally have larger ground sampling distance (GSD) than multispectral (MS) imaging sensors. Data fusion of low spatial-resolution HS and high spatial-resolution MS images enables the production of high spatial-resolution HS data with small spectral distortion [1], [2]. The fused data is useful for accurate classification with fine spatial resolution and thereby enhance applications of HS remote sensing. Several HS-MS fusion methods are proposed using pan-sharpening techniques, stochastic method, and unmixing [1], [2]. A maximum a posteriori (MAP) estimation method was developed to enhance the spatial resolution of HS data using higher spatial-resolution data such as MS and panchromatic images [1]. This approach used a stochastic mixing model (SMM), which estimates the underlying spectral scene characteristics, in order to develop a cost function that optimizes the estimated HS data relative to the observed HS and MS data. An unmixing based fusion method, namely the ’coupled nonnegative matrix factorization’ (CNMF), was recently proposed to enhance spatial resolution of all HS bands. It is based on a linear mixing model (LMM) [2]. This algorithm extracts endmember spectra and high spatial-resolution abundance maps from two images by alternate nonnegative matrix factorization (NMF) [3]. Since CNMF can increase the number of endmembers compared with the MAP/SMM, it has a possibility to deal with spectrally more varied scenes. Many researchers have worked on spectral unmixing with the LMM that assumes that an observed spectrum is a linear combination of several endmember spectra. Although LMM based unmixing methods can obtain physically meaningful results, nonlinearity can appear in spectral mixing model [4], [5]. In recent years, nonlinear unmixing for HS images is receiving growing attention in remote sensing image interpretation. A bilinear mixing model (BMM), which considers second scattering of photons between two distinct materials, has been studied by several groups [6], [7]. A generalized bilinear model (GBM) introduces an effective mean to deal with the underlying assumptions in the BMM [7], [8]. In this work, the authors apply the GBM unmixing to HS-MS data fusion to improve the quality of the fused data introducing semi-nonnegative matrix factorization (Semi-NMF) [9] as a new optimization for the GBM unmixing. 2. GBM VIA SEMI-NMF The BMM considers second-order interactions between different D endmembers as additional terms in the LMM assuming that third or higher order interactions are negligible. In the BMM, the observed L-spectrum of a single pixel z ∈ RL×1 is given by: z = Ea + D−1 ∑ D ∑ i=1 j=i+1 bi,j ei ⊙ ej + n, (1) where E ∈ RL×D is the endmember matrix with the i th column vector, ei ∈ RL×1 representing the i th endmember spectrum, a ∈ RD×1 is the abundance vector, bi,j is the interaction abundance between the i th and j th endmembers, ⊙ is the Hadamard (element-wise product) operation, and n ∈ RL×1 is the additive noise. On the right side, the first term denotes the linear mixing and the second term represents the bilinear mixing. From a physical perspective, the GBM introduces the nonlinear mixing coefficient ci,j as bi,j = ci,j ai aj and assumes the following constraints: ai ≥ 0 ∀i ∈ 1, ..., D and D ∑ ai = 1, i=1 (2) 0 ≤ ci,j ≤ 1 ∀i ∈ 1, ..., D − 1 ∀j ∈ i + 1, ..., D. When the endmembers are known, the GBM unmixing turns to the optimization of ai and ci,j under the constraint of (2). Several optimization methods are proposed in [8]. The GBM method was applied to small images of synthetic and real HS data with three endmembers [7], [8]. When applied to larger images in an unsupervised manner with many endmembers, the optimization process becomes more challenging. The new optimization method based on semi-nonnegative matrix factorization (Semi-NMF) [9] is introduced to speed up the optimization process of a whole image in a matrix form, which can be easily incorporated with CNMF. The observed HS image can be reshaped as a matrix form Z ∈ RL×P with P representing the number of pixels. The BMM for the whole image is given in a matrix form by Z = EA + MB + N, (3) where A ∈ RD×P is the abundance matrix, M ∈ RL×D(D−1)/2 is the bilinear endmember matrix, B ∈ RD(D−1)/2×P is the interaction abundance matrix, and N ∈ RL×P is the noise matrix. The GBM unmixing becomes the ∑ minimization of D ∥Z − EA − MB∥2F , subject to A ⪰ 0 , i=1 Ai,k = 1 ∗ ∗ (∀k ∈ 1, ..., P ), 0 ⪯ B ⪯ A , where A(i,j),k = Ai,k Aj,k . By introducing Z1 = Z − MB and Z2 = Z − EA, (3) is written as follows Z1 = EA + N and Z2 = MB + N. (4) Owing to physical constraints, all components of E, M, A, and B are nonnegative. Therefore, minimization of ∥N∥2F in (4) can be solved by Semi-NMF that factorizes a nonrestricted matrix X into a non-restricted matrix F and a nonnegative matrix GT as X ≈ FGT [9]. Semi-NMF optimization is guaranteed to converge to a local optimum with the alternative update rules. With E given and M calculated from E, the GBM unmixing is solved by the following update rules for A and B: √ T T T T T + − − + (5) AT ←AT .∗ ((ZT 1 E) +A (E E) )./((Z1 E) +A (E E) ) √ T T T T T + − − + (6) BT ←BT .∗ ((ZT 2 M) +B (M M) )./((Z2 M) +B (M M) ) where .∗ and ./ denote elementwise multiplication and division. (C)+ and (C)− are the positive and negative parts of a matrix C defined as C+ = (|C| + C)/2, C− = (|C| − C)/2. The GBM unmixing via Semi-NMF is as follows. First, A is initialized by the fully constrained least square (FCLS) method [10], A∗ is calculated, and B is set as δ × A∗ with small value of δ. Next, A and B are alternately updated by (5) and (6). If any element of B exceeds that of A∗ , it is replaced by that of A∗ . To satisfy the abundance sum-to-one constraint, the method from [10] is adopted. 3. HS-MS DATA FUSION BASED ON GBM The aim of HS and MS data fusion is to estimate unobservable high spatial-resolution HS data (Z ∈ RL×P ) from observed low spatial-resolution HS data (X ∈ RL×Ph ) and high spatial-resolution MS data (Y ∈ RLm ×P ). Lm denotes the number of spectral channels of MS sensor. Ph denotes the number of pixels of HS images. Owing to the trade-off between spectral and spatial resolutions of two sensors, L > Lm and Ph < P are satisfied. The observed two data are assumed to be obtained under the same atmospheric and illumination conditions, and geometrically co-registered with radiometric correction. The HS and MS images can be considered as the degraded versions from the high spatial-resolution HS image in spatial and spectral domains, respectively. Therefore, X and Y are modeled as X = ZS + Ns , (7) Y = RZ + Nr . (8) Here, S ∈ RP ×Ph is the spatial spread transform matrix P ×1 h with each column vector {sl }P representing the l=1 ∈ R transform of the point spread function (PSF) from the MS image to the HS l th∑ pixel value. Each PSF is assumed to Lm Lm ×L be normalized, i.e., is the k=1 skl = 1. R ∈ R spectral response transform matrix with each row vector 1×L m {ri }L representing the transform of the spectral i=1 ∈ R response function (SRF) from the HS sensor to the MS i th band detector. Ns and Nr are the residuals. 3.1. The CNMF algorithm First, we summarize the CNMF method that was proposed for HS-MS data fusion based on the LMM [2]. The CNMF method is composed of alternate NMF unmixing for HS and MS images to extract high spectral-resolution endmember spectra and high spatial-resolution abundance maps. In the LMM, the spectrum at each pixel is assumed to be a linear combination of several endmember spectra, which is the simple version without the bilinear term in (3). Therefore, Z is formulated as Z ≈ EA. The spatially degraded abundance matrix Ah ∈ RD×Ph and the spectrally degraded endmember matrix Em ∈ RLm ×D are defined as Ah = AS and Em = RE, respectively. By substituting Z ≈ EA into (7) and (8), X and Y are approximated as X ≈ EAh and Y ≈ Em A, respectively. GBM for HS data (X) L X ≈ GBM of MS data (Y) P Lm Y D(D-1)/2 D Ph E R ≈ Em Ah M + Bilinear Interpolation A R + Mm Bh Bilinear Interpolation B Fig. 1. Illustration of HS-MS data fusion based on GBM. Owing to the non-negative characteristics of S, R, E, and A, all components of Ah and Em are also non-negative. Z can be estimated by unmixing X and Y to obtain E and A, using alternate NMF unmixing. The CNMF algorithm starts from the NMF unmixing of X to use its spectral advantage. E is initialized by the vertex component analysis (VCA) [11]. Next, X and Y are alternately unmixed by NMF using Ah = AS and Em = RE for the initialization of Ah and Em , respectively. Finally, the high spatial-resolution HS image is estimated by multiplication of E and A. More details concerning CNMF is given in [2]. 3.2. CNMF with GBM The GBM via Semi-NMF can be incorporated into the framework of the CNMF algorithm. Here, the authors define the spatially degraded interaction abundance matrix Bh ∈ RD×Ph and the spectrally degraded bilinear endmember matrix Mm ∈ RLm ×D as Bh = BS and Mm = RM, respectively. By substituting (3) into (7) and (8), X and Y can be approximated as X ≈ EAh + MBh , (9) Y ≈ Em A + Mm B. (10) First, E, Ah , Em , and A are initialized by CNMF. Next, Ah and Bh are estimated by the GBM unmixing of X via Semi-NMF. A and B are initialized by bilinear interpolation of Ah and Bh , respectively. Finally, A and B are optimized by the GBM unmixing of Y, and the fused image (Z) is obtained by Z ≈ EA + MB. The bilinear interpolation is empirically important for the initialization of A and B in the last GBM unmixing of MS image because the GBM unmixing via Semi-NMF converges to a local minimum. Fig. 1 shows the illustration of the proposed method (GBM-CNMF). 4. EXPERIMENTAL STUDY The proposed HS-MS data fusion technique is applied to synthetic datasets generated from the airborne visible infrared imaging spectrometer (AVIRIS) data. The image was taken over Cuprite, Nevada, in 1997 with 224 spectral bands in the 400-2500 nm region. The data initially measured as radiance was converted into reflectance. The 150x150 pixel-size image is selected with 189 bands eliminating noisy bands and the MS and low spatial-resolution HS datasets are generated by down-sampling the original HS data in the spectral and spatial domains, respectively. The spatial-resolution difference and SNR of two sensors are determined considering the specification of Hyperspectral Imager Suite (HISUI) [12], which is composed of HS and MS sensors and will be launched on the Japanese next generation earth observing satellite (Advanced Land Observing Satellite 3 (ALOS-3)) in 2015. The MS data is produced with an uniform SRF corresponding to Landsat TM bands 1-5 and 7. The low spatial-resolution HS data is generated by a Gaussian PSF with full width at half maximum, corresponding to six pixels in the original high spatialresolution HS image, which makes a six times difference in the spatial resolution between two sensors. Therefore, R and S are given as sparse matrices with each row vector corresponding to a uniform SRF and with each column vector representing a Gaussian PSF, respectively. In addition, Gaussian noise was added to the two data assuming that the signal-tonoise ratios (SNR) of HS and MS sensors are 300 and 200, respectively. The performance of HS and MS data fusion was evaluated by comparing the estimated high spatial-resolution HS data with the original data with two criteria: the spatial quality of each spectral band image and the spectral quality of each spectrum for a single pixel. To evaluate the spatial reconstruction quality, the peak signal-to-noise ratio (PSNR) is adopted. To evaluate the spectral reconstruction quality, the authors used the spectral angle error (SAE) between the estimated spectrum and the actual spectrum. Fig. 2 shows the low spatial-resolution HS image and the fused image obtained by GBM-CNMF. The proposed method is evaluated by comparing with CNMF. Fig. 3 shows the average values of PSNR and SAE, and computational time with a varying number of endmembers (D). The performances of two methods converge around D=20 and GBM-CNMF shows better results in the two criteria with all numbers of endmembers. This result proves that the nonlinear mixing effect of multiple reflection is contained in the original HS image and that the GBM unmixing via Semi-NMF works well for both HS and MS images reducing unmixing residual errors. Although there is an obvious quality improvement of GBMCNMF from CNMF, computational time becomes large. For practical use, regarding the computational time, the data fusion process should be processed with partial images using parallel computing. Fig. 4 shows the PSNRs for all bands and the histogram of SAE with SAE maps with D=10. GBMCNMF improves the PSNRs of the fused data in near infrared and short wavelength infrared spectral ranges. It indicates that the GBM unmixing mainly deal with multiple scattering effects in these spectral ranges for this scene. The histogram of Multispectral bands (b) 10 (a) 55 Number of Pixels PSNR (dB) 50 45 40 35 30 25 Fig. 2. Low spatial-resolution HS image (left) and GBMCNMF fused image (right). (a) 44.5 (b) 44 1.05 SAE (deg) PSNR (dB) 42.5 42 41.5 41 GBM-CNMF CNMF 40.5 10 15 20 25 30 35 0.85 0.7 5 10 15 20 25 30 35 40 Number of endmembers 10 5 Computational Time (sec) CNMF 10 1 0 1000 1500 2000 Wavelength (nm) 2500 10 0 0 1 2 3 4 5 SAE (deg) 6. REFERENCES Number of endmembers (c) 500 0.9 0.8 40 10 2 GBM-CNMF CNMF 0.95 0.75 40 39.5 5 4 10 3 GBM-CNMF CNMF 1 43 GBM-CNMF GBM-CNMF CNMF Fig. 4. Comparisons of (a) PSNRs and (b) histogram of SAE with SAE distribution maps between GBM and LMM. 1.1 43.5 4 [1] M. T. Eismann and R. C. Hardie, “Hyperspectral resolution enhancement using high-resolution multispectral imagery with arbitrary response functions,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 3, pp. 455–465, Mar. 2005. [2] N. Yokoya, T. Yairi, and A. Iwasaki, “Coupled nonnegative matrix factorization unmixing for hyperspectral and multispectral data fusion,” IEEE Trans. Geosci. Remote Sens., vol. 50, no. 2, pp. 528–537, Feb. 2012. GBM-CNMF CNMF 10 4 10 3 [3] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, vol. 401, pp. 788-791, Oct. 1999. 10 2 10 1 5 10 15 20 25 30 35 40 Number of endmembers Fig. 3. Effect of the number of endmembers on (a) PSNR, (b) SAE, and (c) computational time. SAE and SAE maps show the spectral quality improvement of GBM-CNMF. The areas where there remain some amount of SAE have a possibility to contain the other nonlinear mixing effect. 5. CONCLUSION In this work, a new HS-MS data fusion method based on the GBM unmixing is proposed. A semi-NMF based optimization method is introduced for the GBM unmixing and incorporated into the HS-MS data fusion algorithm based on linear unmixing (CNMF) to consider the nonlinear mixing caused by multiple scattering. The experiments with synthetic HSMS datasets generated from AVIRIS image show that the proposed HS-MS data fusion method based on the nonlinear unmixing outperforms the CNMF method. The GBM unmixing works well for both HS and MS images with better unmixing residual errors and thereby improves the quality of the HS-MS fused data. Future works include further experiments for various kinds of scenes and the improvement of computational cost using parallel computing. [4] D.A. Roberts, M.O. Smith, and J.B. Adams, “Green vegetation, nonphotosynthetic vegetation, and soils in AVIRIS data,” Remote Sens. Environ., vol. 44, pp. 117–126, 1993. [5] L. Zhang, D. Li, Q. Tong, and L. Zheng, “Study of the spectral mixture model of soil and vegetation in PoYang Lake area, China,” Int. J. Remote Sens., vol. 19, no. 11, pp. 2077–2084, 1998. [6] W. Fan, B. Hu, J. Miller, and M. Li, “Comparative study between a new nonlinear model and common linear model for analysing laboratory simulated-forest hyperspectral data,” Int. J. Remote Sens., vol. 30, no. 11, pp. 2951–2962, Jun. 2009. [7] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, “Nonlinear unmixing of hyperspectral images using a generalized bilinear model,” IEEE Trans. Geosci. Remote Sens., vol. 49, no. 11, pp. 4153–4162, Nov. 2011. [8] A. Halimi, Y. Altmann, N. Dobigeon, and J. Tourneret, “Unmixing hyperspectral images using the generalized bilinear model,” in Proc. IEEE IGARSS, Vancouver, Canada, 2011, pp. 519–522. [9] C. Ding, T. Li, and M. I. Jordan, “Convex and semi-nonnegative matrix factorization,” IEEE Trans. Patten Anal. Mach. Intell., vol. 32, no. 1, pp. 45–55, 2010. [10] D. C. Heinz and C.-I. Chang, “Fully constrained least squares linear spectral mixture analysis method for material quantification in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sens., vol. 39, no. 3, pp. 529–545, Mar. 2001. [11] J. M. Nascimento and J. M. Bioucas-Dias, “Vertex component analysis: A fast algorithm to unmix hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol. 43, no. 4, pp. 898–910, Apr. 2005. [12] N. Ohgi, A. Iwasaki, T. Kawashima, and H. Inada, “Japanese hypermulti spectral mission,” in Proc. IEEE IGARSS, Hawaii, USA, 2010, pp. 3756-3759.