HYPERSPECTRAL AND MULTISPECTRAL DATA FUSIONBASED ON NONLINEAR UNMIXINGNaoto 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) andhigh spatial-resolution multispectral (MS) images based ona linear mixing model (LMM) enables the production ofhigh 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 resultson 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 andhigh spatial-resolution MS images enables the production ofhigh 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 datasuch as MS and panchromatic images [1]. This approach
used a stochastic mixing model (SMM), which estimates theunderlying 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 isreceiving 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 generalizedbilinear model (GBM) introduces an effective mean to deal
with the underlying assumptions in the BMM [7], [8]. In thiswork, 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 isgiven by:z = Ea +

D−1
∑

D
∑

bi,j ei ⊙ ej + n,(1)i=1 j=i+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 andassumes the following constraints:ai ≥ 0 ∀i ∈ 1, ..., D andD∑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 mannerwith many endmembers, the optimization process becomes
more challenging. The new optimization method based onsemi-nonnegative matrix factorization (Semi-NMF) [9] isintroduced 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 amatrix form Z ∈ RL×P with P representing the number ofpixels. 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 theinteraction abundance matrix, and N ∈ RL×P is the noise
matrix. The GBM unmixing becomes the ∑minimization ofD∥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 thealternative update rules. With E given and M calculated fromE, the GBM unmixing is solved by the following update rules
for A and B:
√
T
T
T
T
T
+
−
−
+
AT ←AT .∗ ((ZT(5)
1 E) +A (E E) )./((Z1 E) +A (E E) )
√
T
T
T
T
T
+
−
−
+
BT ←BT .∗ ((ZT(6)
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 amatrix 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∗ withsmall value of δ. Next, A and B are alternately updated by (5)and (6). If any element of B exceeds that of A∗ , it is replacedby 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 highspatial-resolution MS data (Y ∈ RLm ×P ). Lm denotes thenumber of spectral channels of MS sensor. Ph denotes thenumber 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 radiometriccorrection. The HS and MS images can be considered as thedegraded 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 ×1hwith 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 toLm
Lm ×L
is thebe normalized, i.e.,
k=1 skl = 1. R ∈ R
spectral response transform matrix with each row vector
1×Lm{ri }Lrepresenting the transform of the spectral
i=1 ∈ R
response function (SRF) from the HS sensor to the MS i thband 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 isformulated 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
LmY

D(D-1)/2

D

Ph

E

R
≈ Em

Ah

M

+

BilinearInterpolation

A

R
+

Mm

Bh

BilinearInterpolation

B

Fig. 1. Illustration of HS-MS data fusion based on GBM.
Owing to the non-negative characteristics of S, R, E, andA, all components of Ah and Em are also non-negative. Z canbe estimated by unmixing X and Y to obtain E and A, using
alternate NMF unmixing. The CNMF algorithm starts fromthe 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 bymultiplication of E and A. More details concerning CNMF isgiven in [2].
3.2. CNMF with GBM
The GBM via Semi-NMF can be incorporated into theframework 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 canbe 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 interpolationof Ah and Bh , respectively. Finally, A and B are optimizedby 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 theillustration 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 takenover 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 bydown-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 iscomposed of HS and MS sensors and will be launched on theJapanese 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 inthe 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 ofeach 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, theauthors used the spectral angle error (SAE) between the estimated spectrum and the actual spectrum.
Fig. 2 shows the low spatial-resolution HS image and thefused 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 witha 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 ofmultiple reflection is contained in the original HS image andthat the GBM unmixing via Semi-NMF works well for bothHS 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 infraredand 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).
(b) 1.1

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 5Computational 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.

(a) 44.5
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. IEEEIGARSS, 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.