HYPERSPECTRAL, MULTISPECTRAL, AND PANCHROMATIC DATA FUSION BASED ON
COUPLED NON-NEGATIVE MATRIX FACTORIZATION
Naoto Yokoya,1 Takehisa Yairi,2 and Akira Iwasaki2Department of Aeronautics and Astronautics, The University of Tokyo, Japan
Research Center for Advanced Science and Technology, The University of Tokyo, Japan
12ABSTRACT
Coupled non-negative matrix factorization (CNMF) is applied to hyperspectral, multispectral, and panchromatic datafusion. This unmixing based method extracts and fuses hyperspectral endmember spectra and high-spatial-resolution
abundance maps using these three data. An experiment withthe synthetic data simulating ALOS-3 (advanced land observing satellite 3) dataset shows that the CNMF method has apossibility to produce fused data which have both high spatialand spectral resolutions with smaller spectral distortion.
Index Terms— Data fusion, non-negative matrix factorization, spectral unmixing
1. INTRODUCTION
Panchromatic, multispectral, and hyperspectral data analysis
plays an important role in remote sensing. There is always
some engineering trade-off between spatial and spectral resolutions in the optical sensor design. The hyperspectral andmultispectral mission named HISUI (hyperspectral imagersuite) is the Japanese next-generation spaceborne radiometers[1]. The performance of the hyperspectral radiometer is 30m ground sampling distance (GSD), 30 km swath width, and186 spectral channels over 400-2500 nm. The performanceof the multispectral radiometer is 5 m GSD, 90 km swath
width, and 4 spectral channels over 450-900 nm. HISUI will
be launched on advanced land observing satellite 3 (ALOS-3)
of Japan aerospace exploration agency (JAXA) in FY2014.
ALOS-3 also includes the high-spatial-resolution (approximately 1 m GDS) and wide-swath-width (50km) panchromatic sensor with a stereo view. Data fusion of hyperspectral,multispectral, and panchromatic data has a possibility to produce data with both higher spatial and spectral resolutions,which contribute to the accurate identification and classification of the observed materials at fine spatial resolution.
In this work, we propose coupled non-negative matrix
factorization (CNMF) for hyperspectral, multispectral, andpanchromatic data fusion. In recent decades, non-negative
matrix factorization (NMF) has become a useful hyperspectral unmixing method [2]-[4]. Given a non-negative matrix
V, NMF looks for two non-negative matrix factors W andтрс΀к΀норр΀ллйк΀кͿккͿглпͨййЫлйккΒΎΎΎH such that V = WH [5], [6]. We applied NMF to hyperspectral, multispectral, and panchromatic data fusion. Since
it is based on spectral unmixing, the CNMF fused data havelittle spectral distortion while enhancing spatial resolutionof all hyperspectral band images. The CNMF is applied tothe ALOS-3 simulation dataset generated from real airbornehyperspectral data taken over vegetation areas.
2. THE CNMF ALGORITHM
The objective of this work is to estimate unobservablehigh-spatial-resolution hyperspectral data (Z ∈ Rλh ×Lp )from observable low-spatial-resolution hyperspectral (Xh ∈Rλh ×Lh ), multispectral (Xm ∈ Rλm ×Lm ), and panchromatic
data (Xp ∈ R1×Lp ). Here, λh and λm denote the numbers ofspectral channels of hyperspectral and multispectral sensors,respectively. Lh , Lm , and Lp denote the numbers of pixelsof hyperspectral, multispectral, and panchromatic images,respectively.
2.1. Observation Model
The relation between the low- and high-spatial-resolution hyperspectral data can be modeled in matrix form by
Xh = ZSh + Nh ,(1)where Sh ∈ RLp ×Lh is the spatial spread transform maLhtrix with each column vector {sh,k }k=1
representing thetransform from the point spread function of the high-spatialresolution hyperspectral data to that of the k th pixel value
in the low-spatial-resolution hyperspectral data. Nh is theresidual error. Similarly, the multispectral data is modeled as
Xm = Rm ZSm + Nm ,(2)where Rm ∈ Rλm ×λh is the spectral response transformλmmatrix with each row vector {rm,i }i=1representing thetransform from the spectral response function of the hyperspectral data to that of the i th band in the multispectral
data. Sm ∈ RLp ×Lm is the spatial spread transform maLm
trix with each column vector {sm,k }k=1representing thetransform from the point spread function of the high-spatialresolution hyperspectral data to that of the k th pixel value inthe multispectral data. Nm denotes the residual error. The
panchromatic data is also modeled as
Xp = Rp Z + Np ,(3)where Rp ∈ R1×λh is the spectral response transform vector which represents the transform from the spectral response
function of the hyperspectral data to that of the panchromaticdata. Np denotes the residual error. In the simulation of thiswork with synthetic data, Sh , Sm , Rm , and Rp are given.
For real data, Sh and Sm can be determined by image registrations and estimation of point spread functions. Rm andRp can be also derived from radiometric calibration to obtain
spectral response functions.
2.2. Coupled NMF Unmixing
With the linear spectral mixture assumption, Xh , Xm and Xp
are expressed as follows:
Xh = Wh Hh + Eh ,(4)

Xm = Wm Hm + Em ,(5)

Xp = Wp Hp + Ep .

(6)

Here, Wh ∈ Rλh ×D , Hh ∈ RD×Lh , and Eh ∈ Rλh ×Lh
are the endmember, abundance, and residual matrices of the
low-spatial-resolution hyperspectral data, respectively. D isthe number of endmember. Wm ∈ Rλm ×D , Hm ∈ RD×Lm ,and Em ∈ Rλm ×Lm are those of the multispectral data. Inthe same way, Wp ∈ R1×D , Hp ∈ RD×Lp , and Ep ∈ R1×Lp
are those of the panchromatic data. It is physically reasonable to assume that the high-spatial-resolution hyperspectral
data contains the same endmember spectra as the low-spatialresolution hyperspectral data and the same abundance maps
as the panchromatic data. Therefore, Z can be approximatedasZ ≈ Wh Hp .
(7)
This is the key idea of unmixing based data fusion. CNMF iscomposed of two steps based on NMF, i.e., 1) the unmixingsof Xh and Xm , and 2) the estimation of Hp using Wh and Hm .
Finally, by combining Wh and Hp , the high-spatial-resolution
hyperspectral data can be obtained. Other unmixing and endmember detection algorithms are also applicable for this data
fusion. But, to minimize the residual errors in the spectral
mixture models taking account of sensor properties, NMF isstraightforward to formulate and easy to implement.
2.2.1. Unmixings of hyperspectral and multispectral data
First, Xh and Xm are alternately unmixed by NMF to estimate Wh and Hm . NMF spectral unmixing is commonly performed to minimize the squared Frobenius norm of the residual matrix in the linear spectral mixture model expressed as∥Eh ∥2F and ∥Em ∥2F for (4) and (5), respectively. Lee and Seung proposed a multiplicative update rule that is guaranteed
to converge to a locally optimal matrix factorization under
the non-negative constraints of two factorized matrices. We
use the multiplicative update rule for NMF unmixings of Xhand Xm . From (1), (2), (4), (5), and (7), the endmember andabundance matrices are related asS∗hHh ≈ Hm S∗h ,(8)Wm ≈ Rm Wh .

(9)is the spatial spread transform matrix with"L h!
each column vector s∗h,krepresenting the transform
k=1
from the point spread function of the multispectral data tothat of the k th pixel value in the hyperspectral data, which
can be obtained from Sh and Sm . Xh and Xm are alternately
unmixed to estimate Wh and Hm under constraints of (8) and(9).
The CNMF algorithm starts from NMF of Xh to use its
spectral advantage. First, with the number of endmember D
set to a certain number, Wh is initialized by vertex component
analysis (VCA) [7], which is one of the most advanced geometry based endmember extraction methods with pure pixel assumption. Hh is initialized with a constant value 1/D. Next,
Wh and Hh are optimized by the multiplicative update rule.
However, as the initialization phase, only Hh is updated until convergence with Wh fixed. In the subsequent rounds, thevalue of Hh that is initialized by (8) is used and Wh is updated with Hh fixed until convergence to inherit the reliable
abundance information obtained from NMF of Xm .
The alternate step is NMF of Xm . Wm is initialized by(9) and Hm is set as a constant value 1/D. As the initialization phase, Hm is updated until convergence with Wm fixed
to inherit the reliable endmember spectra obtained from NMFof Xh . Next, Wm and Hm are optimized by the multiplicative update rule. We refer to the alternate NMF unmixings as
the outer loops and the update processes in each NMF unmxing as the inner loops. As a convergence condition, we use
the condition that the change ratio of cost function achieves
a value below a given threshold. Considering practical use, amaximum number of iterations is set for each loop. Two parameters are set for the outer and inner loops, respectively, asdifferent values.
∈ R

Lm ×Lh2.2.2. Estimation of high-spatial-resolution abundance matrix from panchromatic data
In the last step, the high-spatial-resolution abundance matrix
is estimated from the panchromatic data by NMF. Since thepanchromatic data have only one band, the initial value of Wp
and Hp is significantly important for the accurate NMF of Xp .
From (3), (6), and (7), Wp is approximated asWp ≈ Rp Wh .

(10)

We estimate Wp by (10) using Wh obtained in the previousphase. Hp are initialized by the spatially up-sampled matrixof Hm using bilinear interpolation. Only Hp is optimized because the panchromatic data have no spectral information. To
inherit the abundance information of Hm , we minimize thefollowing cost function: ∥Ep ∥2F + α∥Hp − H̃m ∥2F , whereH̃m ∈ RD×Lp denotes bilinearly up-sampled Hm and α isa penalty parameter. The multiplicative update rule of Hp isgiven by472nm

765nm

Zoom (765nm)(a)90m(b)Hp ← Hp .∗(Wp T Xp −α(Hp − H̃m ))./(Wp T Wp Hp ), (11)where (·)T denotes the transposition of the matrix and .∗ and./ denote element-wise multiplication and division, respectively. Hp is optimized by (11) until convergence. Finally, thehigh-spatial-resolution hyperspectral data are obtained by themultiplication of Wh and Hp .
3. DATA AND METRIC(d)The proposed data fusion technique is applied to the ALOS-3
simulation dataset generated from real airborne hyperspectral
data. The image was taken over a vegetation area in Japanby the compact airborne spectrographic imager (CASI) sensor
with 68 spectral bands in the 400-1060 nm region in 2008. We
selected a 480×480 pixel size image with 47 bands in the 450900 nm region and generated the panchromatic, multispectral,and low-resolution hyperspectral data by down-sampling theoriginal hyperspectral data in the spectral, both spectral andspatial, and spatial domains, respectively.
The panchromatic data was produced with a uniform
spectral response function in the 450-900 nm region. The
multispectral data was produced by averaging 5×5 pixel
blocks of the original high-resolution hyperspectral data withuniform spectral response functions corresponding to multispectral bands 1-4 of HISUI, which cover the 450-520,520-600, 630-690, and 760-900 nm regions, respectively [1].
The low-resolution hyperspectral data was generated by averaging 30×30 pixel blocks of the original high-resolution
hyperspectral data. In this simulation, we assume completelyidentical atmospheric and illumination conditions between
the three sensors.
We evaluate the spatial and spectral qualities of the fused
data compared with the original data using the peak signalto-noise ratio (PSNR) and the spectral angle error (SAE), respectively. The PSNR of the i th band is defined byPSNRi = 10 · log10#max2iLpΣk=1(Z − Wh Hp )2i,k /Lp(c)$, (12)where maxi is the maximum pixel value in the i th band image. The SAE of each pixel is defined as an angle in λh spectral dimensions.

(e)Fig. 1. 472 nm, 765 nm, and 765 nm zoomed images of (a)low-spatial-resolution hyperspectral, (b) original, (c) HM, (d)HP, and (e) HMP data.
4. RESULTS AND DISCUSSION
We compare the quality of three fused data produced byCNMF, i.e., the hyperspectral and multispectral (HM), hyperspectral and panchromatic (HP), and hyperspectral, multispectral, and panchromatic (HMP) fused data. HM is obtained in the first phase of CNMF. HP can be produced in thesecond phase of CNMF substituting Hh for Hm in the initialization of Hp . Fig. 1 shows the comparison between thelow-spatial-resolution hyperspectral (HS), the original, andthe three fused data in the 472 nm and 765 nm regions. In thetwo band images, HMP is most similar to the original image
regarding high spatial resolution. In the 765 nm band image,HP seems to be comparable with HMP. However, in the 472
nm band image, the reflectance values are reversed in many
areas. This indicates that it is difficult to estimate abundancemaps with high spatial resolution (1 m) from those with low
spatial resolution (30 m).

Multispectral
Panchromatic
40(a)

36

1034
PSNR

10

Number of Pixels

HMHP
HMP

38

1032

1030
28

1026
24
22
400

10500

600
700
800
Wavelength (nm)900

5HMHPHMP432

100

102030SAE40 (deg)(b)

Fig. 2. PSNR of HM, HP, and HMP fused data.
Fig. 2 shows the PSNR of the three fused data. HMP outperforms HM and HP in all spectral channels. This provesthat both multispectral and panchromatic data contribute tobetter estimation of the high-spatial-resolution hyperspectral
data. Fig. 3 shows the histograms of SAE and the SAE distribution maps. HP shows large spectral distortion due to thelack of spectral information in the panchromatic data. The result that HM has less spectral distortion than HP means thatthe multispectral data are useful for resolution enhancementof the low-spatial-resolution hyperspectral data while keeping
the spectral distortions small. HMP is thought to be comparable with HM because the panchromatic data do not contribute to the improvement of spectral quality. As shown in
the SAE distribution map of HMP, SAE values are small inalmost all vegetation areas and spectral distortions occur in
non-vegetation area. The reason is as follows. Endmember
matrix (Wh ) mainly contains vegetation spectra because thelow-spatial-resolution hyperspectral data has a large coverageof vegetation area whose spectra is highly mixed. Therefore,
it is difficult to linearly combine the endmembers to represent
the spectra of non-vegetation materials. If we can estimate
more accurate endmember spectra in the first step of CNMF,
the mitigation of spectral distortions may be possible.
5. CONCLUSION
We proposed the CNMF algorithm for the fusion of panchromatic, multispectral, and hyperspectral data. First, the multispectral and hyperspectral data are alternately unmixed into
the endmember and abundance matrices by NMF with therelation between sensor properties used for the constraints.
Next, the high-spatial-resolution abundance maps are estimated from the panchromatic data using the hyperspectral
endmember and multispectral abundance matrices. Finally,
the high-spatial-resolution hyperspectral data is obtainedby multiplication of the hyperspectral endmember and highspatial-resolution abundance matrices. Owing to an algorithm
based on unmixing, the CNMF method has a possibility toHMHP

HMPFig. 3. (a) Histograms of SAE and (b) SAE distribution maps
extract and fuse all spectral and spatial information containedin the three data. It is applicable for the high-order data
product of the ALOS-3 dataset.
6. ACKNOWLEDGMENT
The authors would like to acknowledge ERSDAC who provided CASI data.
7. REFERENCES[1] N. Ohgi, A. Iwasaki, T. Kawashima, and H. Inada, “Japanese
hyper-multi spectral mission,” in Proc. IEEE IGARSS, Hawaii,
USA, 2010, pp. 3756-3759.
[2] L. Miao and H. Qi, “Endmember extraction from highly mixed
data using minimum volume constrained nonnegative matrix
factorization,” IEEE Trans. Geosci. Remote Sens., vol. 45, no.
3, pp. 765-777, Mar. 2007.
[3] S. Jia and Y. Qian, “Constrained nonnegative matrix factorization for hyperspectral unmixing,” IEEE Trans. Geosci. RemoteSens., vol. 47, no. 1, pp. 161-173, Jan. 2009.
[4] A. Huck, M. Guillaume, and J. Blanc-Talon, “Minimum dispersion constrained nonnegative matrix factorization to unmix hyperspectral data,” IEEE Trans. Geosci. Remote Sens., vol. 48,
no. 6, pp. 2590-2602, Jun. 2010.
[5] D. D. Lee and H. S. Seung, “Learning the parts of objects bynonnegative matrix factorization,” Nature, vol. 401, pp. 788791, Oct. 1999.
[6] D. D. Lee and H. S. Seung, “Algorithms for non-negative matrix
factorization,” in Proc. Conf. Adv. Neural Inf. Process. Syst.,
2001, vol. 13, pp. 556-562.
[7] J. M. P. Nascimento and J. M. B. Dias, “Vertex component
analysis: A fast algorithm to unmix hyperspectral data,” IEEETrans. Geosci. Remote Sens., vol. 43, no. 4, pp. 898-910, Apr.
2005.