COUPLED NON-NEGATIVE MATRIX FACTORIZATION (CNMF) FOR HYPERSPECTRALAND MULTISPECTRAL DATA FUSION: APPLICATION TO PASTURE CLASSIFICATION
Naoto Yokoya,1 Takehisa Yairi,2 and Akira Iwasaki212Department of Aeronautics and Astronautics, The University of Tokyo, Japan
Research Center for Advanced Science and Technology, The University of Tokyo, Japan
ABSTRACT

Coupled non-negative matrix factorization (CNMF) is introduced for hyperspectral and multispectral data fusion. The
CNMF fused data have little spectral distortion while enhancing spatial resolution of all hyperspectral band images owingto its unmixing based algorithm. CNMF is applied to the synthetic dataset generated from real airborne hyperspectral datataken over pasture area. The spectral quality of fused data isevaluated by the classification accuracy of pasture types. The
experiment result shows that CNMF enables accurate identification and classification of observed materials at fine spatialresolution.
Index Terms— Non-negative matrix factorization, unmixing, data fusion, pasture classification
1. INTRODUCTION
The spatial resolution of hyperspectral sensors is often lower
than that of multispectral sensors. The hyperspectral and multispectral mission named HISUI (Hyperspectral Imager Suite)is the Japanese next-generation spaceborne radiometer [1].
HISUI is composed of a hyperspectral radiometer with 30 m
ground sampling distance (GSD) and 186 spectral channels
over 400-2500 nm, and a multispectral radiometer with 5 m
GSD and 4 spectral channels over 450-900 nm. Data fusion ofhyperspectral and multispectral data has a possibility of producing fused data with both higher spatial and spectral resolutions, which contribute to the accurate identification andclassification of observed materials at fine spatial resolution.
In low-spatial-resolution hyperspectral data, an observed
spectrum at each pixel is a mixture of several endmember
spectra, and is commonly called as “mixel” (mixed pixel).
Unmixing is a technique to decompose the observed spectrum
into endmember spectra and their abundance fractions. In therecent decade, non-negative matrix factorization (NMF) [2],[3], has become one of the most useful hyperspectral unmixing techniques based on a linear spectral mixture model [4][6]. Given a non-negative matrix V, NMF looks for two nonnegative matrix factors W and H such that V = WH. When
NMF is applied to hyperspectral unmixing, V, W, and H correspond to observed data, endmember spectra, and abundance

fractions, respectively.
We proposed coupled non-negative matrix factorization(CNMF) for the fusion of spectroradiometer datasets [7].
In the case of hyperspectral and multispectral data fusion,
CNMF unmixes the low-spatial-resolution hyperspectral data
and the high-spatial-resolution multispectral data alternately
using the relation between sensor properties for the NMF
initializations. By combining the hyperspectral endmember matrix and the high-spatial-resolution abundance matrix,
high-spatial-resolution hyperspectral data that inherit bothadvantages of the two data can be synthesized. The CNMF
method is physically straightforward and easy to implement,
owing to its simple update rules. Since CNMF is based on
spectral unmixing, the CNMF fused data have little spectral
distortion while enhancing spatial resolution of all hyperspectral band images.
In this work, we introduce the CNMF algorithm for hyperspectral and multispectral data fusion and investigate theutility of the fused data for pasture classification. CNMF isapplied to the synthetic dataset generated from real airborne
hyperspectral data taken over pasture area. After numericalevaluation of the spatial and spectral qualities of the fused
data, its utility for pasture classification is demonstrated.

2. COUPLED NON-NEGATIVE MATRIXFACTORIZATION (CNMF)
The aim of hyperspectral and multispectral data fusion
is to estimate unobservable high-spatial-resolution hyperspectral data (Z ∈ Rλh ×Lm ) from observable low-spatialresolution hyperspectral (X ∈ Rλh ×Lh ) and multispectral(Y ∈ Rλm ×Lm ) data, which are spatially co-registered. Here,λh and λm denote the numbers of spectral channels of hyperspectral and multispectral sensors, respectively. Lh and Lm
denote the numbers of pixels of hyperspectral and multispectral data, respectively. λh > λm and Lh < Lm are satisfied
due to the trade-off between spatial and spectral resolutionsof two sensors.

2.1. Observation Model
The relation between the low- and high-spatial-resolution hyperspectral data can be expressed in matrix form byX = ZS + Ns ,This is the key idea of CNMF data fusion. From (1)-(4) and(5), the endmember and abundance matrices are related as(1)

where S ∈ RLm ×Lh is the spatial transform matrix with each
Lhcolumn vector {sk }k=1
representing the transform from thepoint spread function of the high-spatial-resolution hyperspectral data to that of k th pixel value in the low-spatialresolution hyperspectral data. Ns is the residual error. Similarly, the multispectral data is related to the high-spatialresolution hyperspectral data byHh ≈ Hm S,(6)

Wm ≈ RWh .

(7)

X = Wh Hh + Eh ,(3)

CNMF alternately unmixes X and Y to estimate Wh and Hmunder constraints of (6) and (7).
The CNMF algorithm starts from NMF of X to use its
spectral advantage. As the initialization phase, with the number of endmembers D set to a certain value, Wh is initialized by vertex component analysis (VCA), which is one of themost advanced geometry based endmember extraction methods with pure pixel assumption. Hh is set to a constant value
1/D and updated by MUR with Wh fixed. As the optimization phase, Wh and Hh are updated by MUR until convergence. In the subsequent rounds, only the initialization phaseis different, i.e., the value of Hh that is initialized by (6) isused and Wh is updated by MUR with Hh fixed to inherit thereliable abundance information obtained from NMF of Y.
The alternate step of CNMF is NMF of Y. As the initialization phase, Wm is initialized by (7). Hm is set to a constantvalue 1/D and updated by MUR with Wm fixed to inherit thereliable endmember spectra obtained from NMF of X. As theoptimization phase, Wm and Hm are updated by MUR untilconvergence.
We refer to the alternate NMF unmixings as the outerloops and the update processes in each NMF nimxing as theinner loops. As a convergence condition, we use the condition
that the change ratio of cost function achieves a value below agiven threshold. Considering practical use, a maximum number of iterations is set for each loop. Two parameters are set
for the outer and inner loops, respectively, as different values.

Y = Wm Hm + Em .

(4)

3. EXPERIMENTY = RZ + Nr ,(2)

where R ∈ Rλm ×λh is the spectral transform matrix withλmeach row vector {ri }i=1
representing the transform from thespectral response function of the hyperspectral data to the i thband in the multispectral data. Nr denotes the residual error.
In the simulation of this work with synthetic data, S and R aregiven.
When applied to real data, S is determined by image registrations and estimation of point spread functions and R is derived from radiometric calibration to obtain spectral response
functions.
2.2. Coupled NMF Unmixing
With the linear spectral mixture assumption, X and Y are expressed as follows:

Here, Wh ∈ Rλh ×D , Hh ∈ RD×Lh , and Eh ∈ Rλh ×Lh are
the endmember, abundance, and residual matrices of the lowspatial-resolution hyperspectral data, respectively. Wm ∈Rλm ×D , Hm ∈ RD×Lm , and Em ∈ Rλm ×Lm are those of
the multispectral data. NMF spectral unmixing is commonly
performed to minimize the squared Frobenius norm of theresidual matrix in the linear spectral mixture model expressedas ∥Eh ∥2F and ∥Em ∥2F for (3) and (4), respectively. Lee andSeung proposed a multiplicative update rule (MUR) that is
guaranteed to converge to a locally optima under the nonnegativity constraints of factorized two matrices. We used
MUR for NMF unmixings of X and Y. 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 multispectral data. Therefore, Z can be approximatedasZ ≈ Wh Hm .
(5)

3.1. Dataset
CNMF is applied to the synthetic dataset generated from compact airborne spectrographic imager 3 (CASI-3) data takenover pasture area in Hokkaido, Japan, on June 19, 2009. The
original data have 34 spectral channels over 410-1070 nmwith 1 m GSD. The data, originally measured as radiance,
was converted into reflectance. By down-sampling the original high-spatial-resolution hyperspectral data in the spectral
and spatial domains, we generate the synthetic dataset simulating the characteristics of HISUI except for spectral resolution of hyperspectral sensor. We ideally assume rectangle
functions for the point spread and spectral response functionsof two sensors. The multispectral data are produced by averaging 5×5 pixel blocks of the original data with uniform
spectral response functions corresponding to the multispectral bands 1-4 of the HISUI, which cover the 450-520, 520600, 630-690, and 760-900 nm regions, respectively. The

+

1A
0.7

OG
RCG
KB
SVG
ER

0.4
0.2

100m0OG

RCG

KB

SVG

ER500600

B
Reference

Reflectance

0.60.6+

Reflectance

0.8

0.50.4

Reference
Hyperspectral
Multispectral
CNMF

0.30.2700 800 900 1000
Wavelength (nm)0.1

Fig. 1. (a) RGB color image of test area with ground truthclassification samples plotted and (b) spectra of five plants.

A

0400 500 600 700 800 900 1000 1100
Wavelength (nm)Hyperspectral

0.8

Table 1. SAD (deg) between five plants.
OG RCG
KB
SVG
ER
OG
0
2.688 1.008 2.892 2.324
RCG
03.266 5.518 0.906
KB
02.516 2.956
SVG
05.088RE
0Multispectral

Reflectance

0.7
0.6
0.5
0.4
0.3
0.20.1

CNMF

B0400 500 600 700 800 900 1000 1100
Wavelength (nm)Fig. 2. (a) Endmember spectra and continuum removed spectra in (b) 545-760 nm and (c) 900-1035 nm.
low-spatial-resolution hyperspectral data are generated by averaging 30×30 pixel blocks of the original data. In the sameway, we made the hyperspectral data with 5 m GSD as thereference data that we aim to produce by CNMF.
3.2. Evaluation of fused data
The qualities of the CNMF fused data are evaluated by comparing it with the reference data. Peak-signal-to-noise ratio(PSNR) is used to evaluate the spatial reconstructivity of eachband image. Higher PSNR indicates better spatial reconstructivity. Brovey transform, which is one of the most popularpan-sharpen algorithms, is used as the comparison method.
This method is applicable only for spectral regions coveredby the multispectral sensor. For each pixel spectrum, spectral
angle distance (SAD) is used to evaluate the spectral distortion. Smaller SAD indicates less spectral distortion.
3.3. Classification
The utility of the CNMF fused data for classification is evaluated by pasture classification. Observed pasture is mainlycomposed of five plants such as orchard grass (OG), reed canary grass (RCG), Kentucky bluegrass (KB), sweet vernalgrass (SVG), and Elytrigia repens (ER). Ground truth data,i.e., class labels and spectra, were obtained at 45 quadrateareas with 3×3 m size. Analytical Spectral Device (ASD)
FieldSpec Pro spectroradiometer was used for field spectrum
measurement. Fig. 1 shows the RGB color image of the testarea with ground truth classification samples superposed andthe spectra of the five plants. SAD between these spectra aresummarized in Table 1. Pairs of OG and KB, and RCG andER are spectrally similar with about one degree of SAD. Since

the CASI data have a geometrical error of about 2 m, the classification result is considered as correct if the major plant in
a 3×3 pixel square with its center set within 2 m of the sampling point matches the ground truth data. We use the spectral angle mapper (SAM) as a classifier. Especially for thehyperspectral data, the continuum removal method with the545-760 nm and 900-1035 nm regions is adopted [8]. Classification accuracy of the CNMF fused data is compared withthose of the hyperspectral, multispectral, and reference data.
4. RESULTS AND DISCUSSION
Left images of Fig. 2 show RGB image of the CNMF fused
data with those of the reference, hyperspectral, and multispectral data. Right graphs show the comparison of spectra at twopoints. The CNMF fused data is close to the reference data
both in terms of spatial and spectral domains. Fig. 3 (a) shows
the PSNR of the fused data obtained by CNMF and Brovey.
CNMF outperforms Brovey in all spectral bands owing to its
unmixing based optimization algorithm. In the spectral bands
not covered by multispectral data, the PSNR values of theCNMF fused data are relatively low. This indicates that resolution enhancement is limited without high-spatial-resolution
information. Fig. 3 (b) and (c) show the histogram of SAD
and the SAD distribution map, respectively. In a lot of vegetation areas, SAD is around one degree, which is comparable with the minimum SAD between five plants. Streams
and trails show larger SAD. Since there are no pure pixels ofthese objects in the low-spatial-resolution hyperspectral data,
it is difficult to retrieve their endmember spectra by NMF un-Multispectral bands
CNMFBrovey

PSNR55(b)504540(c)(d)Average: 1.35 (deg)Number of Pixels

6035(a)10 56510 410310 2(a)(b)30400 500 600 700 800 900 1000 1100
Wavelength (nm)

10 10246SAD (deg)81010842(c)SAD (deg)

60OG

RCG

KB

SVG

ER

Fig. 4. Classification of pasture types generated from four
data, i.e., the (a) hyperspectral, (b) multispectral, (c) CNMF,and (d) reference data.

Fig. 3. (a) PSNR, (b) histogram of SAD, and (c) SAD distribution map.
mixing. Spectra of water and soil can not be expressed aslinear combinations of vegetation spectra, which results inlarger SAD. To produce fused data with little spectral distortion, endmember extraction is the most important process.
Fig. 4 shows the classification result of the observed area
and Table 2 shows classification accuracies for four data, i.e.,the hyperspectral, multispectral, CNMF, and reference data.
Low classification accuracies of the hyperspectral and multispectral data suggest that both spatial and spectral resolutions have significant effects on accurate classification. TheCNMF fused data shows comparable classification accuracy
with the reference data, outperforming the hyperspectral andmultispectral data. This indicates that the CNMF method has
a potential to produce high-spatial-resolution hyperspectral
data with little spectral distortion that enable accurate identification and classification of observed materials at fine spatialresolution.
5. CONCLUSION
In this work, we introduced CNMF for hyperspectral and multispectral data fusion and demonstrated its utility for accurateclassification. Two data with trade-off between spatial andspectral resolutions are reshaped as matrices and unmixed
into endmember and abundance matrices by NMF. By multiplying the endmember matrix derived from hyperspectral
data unmixing and the high-spatial-resolution abundance matrix derived from multispectral data unmixing, high-spatialresolution hyperspectral data can be obtained. CNMF is applied to the synthetic dataset generated from CASI data takenover pasture areas. Classification experiment of pasture types
proved that CNMF fused data enable accurate identificationand classification of observed materials at fine spatial resolution. This technique is useful for the generation of higherorder data products in the HISUI dataset.

Table 2. Classification accuracies (%) for four data.
Hyperspectral Multispectral CNMF Reference
46.7
51.1
64.4
66.7

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 hypermulti spectral mission,” in Proc. IEEE IGARSS, Jul. 2010, pp. 37563759.
[2] D. D. Lee and H. S. Seung, “Learning the parts of objects by nonnegative matrix factorization,” Nature, vol. 401, pp. 788-791, Oct. 1999.
[3] 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.
[4] L. Miao and H. Qi, “Endmember extraction from highly mixed datausing minimum volume constrained nonnegative matrix factorization,” IEEE Trans. Geosci. Remote Sens., vol. 45, no. 3, pp. 765-777,Mar. 2007.
[5] S. Jia and Y. Qian, “Constrained nonnegative matrix factorization for
hyperspectral unmixing,” IEEE Trans. Geosci. Remote Sens., vol. 47,no. 1, pp. 161-173, Jan. 2009.
[6] 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. 25902602, Jun. 2010.
[7] N. Yokoya, T. Yairi, and A. Iwasaki, “Hyperspectral, multispectral,and panchromatic data fusion based on coupled non-negative matrix
factorization,” in Proc. IEEE WHISPERS, Lisbon, Jun. 2011.
[8] ERSDAC, “Annual report in the research and development of usage
technology for the next generation advanced earth observing satellite,”(in Japanese), pp. 47-51, Mar. 2010.