1

Hyperspectral Image Classification with Markov
Random Fields and a Convolutional Neural Network

arXiv:1705.00727v2 [cs.CV] 12 Nov 2017

Xiangyong Cao† , Feng Zhou, Member, IEEE, Lin Xu, Deyu Meng, Member, IEEE, Zongben Xu,
and John Paisley Member, IEEE

Abstract—This paper presents a new supervised classification
algorithm for remotely sensed hyperspectral image (HSI) which
integrates spectral and spatial information in a unified Bayesian
framework. First, we formulate the HSI classification problem
from a Bayesian perspective. Then, we adopt a convolutional
neural network (CNN) to learn the posterior class distributions
using a patch-wise training strategy to better use the spatial
information. Next, spatial information is further considered by
placing a spatial smoothness prior on the labels. Finally, we
iteratively update the CNN parameters using stochastic gradient
decent (SGD) and update the class labels of all pixel vectors
using α-expansion min-cut-based algorithm. Compared with
other state-of-the-art methods, the proposed classification method
achieves better performance on one synthetic dataset and two
benchmark HSI datasets in a number of experimental settings.

I. I NTRODUCTION
Hyperspectral remote sensors capture digital images in hundreds of continuous narrow spectral bands to produce a highdimensional hyperspectral image (HSI). Since HSI provides
detailed information on spectral and spatial distributions of
distinct materials [49], it has been used for many applications, such as land-use mapping [10], [34], [48], land-cover
mapping [34], forest inventory [43], and urban-area monitoring [53]. All these applications require the material class label
of each hyperspectral pixel vector and thus HSI classification
has been an active research topic in the field of remote
sensing. The aim of HSI classification is to categorize each
hyperspectral pixel vector into a discrete set of meaningful
classes according to the image contents.
In the last few decades, many methods have been proposed
for HSI classification. These methods can be roughly divided
into two categories: spectral based methods and spectralspatial based methods. We first briefly review these HSI
classification approaches and then discuss the contributions
of our proposed method.
Xiangyong Cao, Deyu Meng, and Zongben Xu are with the School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an 710049, China (caoxiangyong45@gmail.com, dymeng@mail.xjtu.edu.cn, zbxu@mail.xjtu.edu.cn).
Feng Zhou is with the National Laboratory of Radar Signal Processing,
Xidian University, Xi’an, China (fzhou@mail.xidian.edu.cn).
Lin Xu is with the NYU Multimedia and Visual Computing Lab, New York
University Abu Dhabi, UAE (xulinshadow@gmail.com).
John Paisley is with the Department of Electrical Engineering &
Data Science Institute, Columbia University, New York, NY, USA (jpaisley@columbia.edu).
This research was supported in part by National Basic Research Program
(973 Program) of China under Grant No. 2013CB329404 the National
Natural Science Foundation of China under Grants 91330204, 11131006 and
61373114.
†Xiangyong Cao performed this work at Columbia University on a grant
from the China Scholarship Council.

A. Related work: spectral vs spectral-spatial based methods
Many classical HSI classification approaches are only based
on spectral information [2], [40], [62]. Among these methods,
pure spectral classification without any band reduction has
often been proposed in the literature [6], [11], [18], [44].
Besides, some other pure spectral methods which extract spectral feature first using some feature extraction methods, such
as principal component analysis [40], independent component
analysis [62] and linear discriminant analysis [2], have also
been proposed. However, these approaches only consider the
spectral information and ignore the correlations among distinct
pixels in the image, which tends to decrease their classification
performance relative to those which consider both. In this
paper we instead focus on designing a spectral-spatial based
classification method.
Spectral-spatial based methods can help improve the classification performance since they incorporate additional spatial
information from the HSI. It has been often observed that
spatial information is often as crucial as spectral information
in the HSI classification task [23], [49], [50], [59]. Therefore,
spectral-spatial methods have been proposed that additionally
consider spatial correlation information. One approach is to
extract the spatial dependence in advance using a spatial feature extraction method before learning a classifier. The patchbased feature extraction method is a representative example
of this approach [15], [16], [55], [57]. Here, features are
extracted from groups of neighboring pixels using a subspace
learning technique such as low-rank matrix factorization [64],
dictionary learning [15], [55] or subspace clustering [32].
Compared with the original spectral vector, the features extracted by patch-based methods have higher spatial smoothness
with some reduced noise [12].
Another popular spectral-spatial approach to incorporate
spatial information is using a Markov random field (MRF)
to post-process the classification map. MRF is an undirected
graphical model [4] that has been applied in a variety of fields
from physics to computer vision and machine learning. In
particular, they have been widely used for image processing
tasks such as image registration [70], image restoration [5],
image compression [51] and image segmentation [13]. In
the image segmentation task, MRFs encourage neighboring
pixels to have the same class label [39]. This has been
shown to greatly improve the classification accuracy in HSI
classification task [12], [38], [61].
Another family of spatial-spectral methods is based on
spatial regularization (anisotropic smoothing) prior to spec-

2

tral classification of the regularized image [22], [45]. Aside
from these methods, some other spatial-spectral methods have
also been proposed, such as 3-dimensional discrete wavelet
transform [12], [50], 3-dimensional Gabor wavelets [54], morphological profiles [3], attribute profiles [20] and manifold
learning methods [19], [42].
Although the previous approaches perform well, a drawback
to many of these approaches is that the extracted features are
hand-crafted. Specifically, these approaches extract the features of the HSI by pre-specified strategies manually designed
directly on data without using the label information and not via
an “end-to-end” manner such as the deep learning approaches
introduced here. Thus they highly depend on prior knowledge of the specific domain and are often sub-optimal [68].
In the last few years, deep learning has been a powerful
machine learning technique for learning data-dependent and
hierarchical feature representations from raw data [25]. Such
methods have been widely applied to image processing and
computer vision problems, such as image classification [37],
image segmentation [13], action recognition [31] and object
detection [58].
Recently, a few deep learning methods have been introduced
for the HSI classification task. For example, unsupervised
feature learning methods such as stacked autoencoders [14]
and deep belief networks [17] have been proposed. Although
these two unsupervised learning models can extract deep
hierarchical feature, the 3-dimensional (3D) patch must be first
flattened into 1-dimensional (1D) vectors in order to satisfy
the input requirement, thus losing some spatial information.
A supervised autoencoder [41] method that used label information during feature learning was also proposed. Similar
methods have been proposed based on the convolutional neural
network (CNN), for example, using a five-layer CNN trained
on a 1D spectral vector without using spatial information [29],
using a modified spectral-spatial deep CNN (SS-DCNN) [67]
with 3D patch input, or deep CNNs with spatial pyramid
pooling (SPP) [27] (SPP-DCNN) [66]. Other methods based
on CNNs [47], [68] and recurrent neural networks (RNNs) [46]
have also been proposed.
B. Contributions of our approach
As mentioned, compared with the traditional spectral-spatial
HSI classification methods, deep learning can directly learn
data-dependent and hierarchical feature representation from
raw data. Although all the above deep learning methods
obtain good performance, none of them formulate the HSI
classification task into a Bayesian framework, where deep
learning and MRF are considered simultaneously. Therefore,
in this paper we propose a new supervised HSI classification
algorithm in a Bayesian framework based on deep learning
and MRF. Our contributions are threefold:
1) We formulate the HSI classification problem from a
Bayesian perspective and solve this problem by introducing
an intermediate variable and other simplifications.
2) We propose a new model which combines a CNN
with a smooth MRF prior. Specifically, the CNN is used
to extract spectral-spatial features from 3D patches and the

smooth MRF prior is placed on the labels to further exploit
spatial information. Optimizing the final model can be done
by iteratively updating CNN parameters and class labels of
all the pixel vectors. In this way, we integrate the MRF with
the CNN. To our knowledge, this is the first approach that
integrates the MRF with the CNN for the HSI classification
problem.
3) Our experimental results on one synthetic HSI dataset
and two real HSI datasets demonstrate that the proposed
method outperforms other state-of-the-art methods for the HSI
classification problem, in particular the three deep learning
methods SS-DCNN [67], SPP-DCNN [66] and DC-CNN [68].
The rest of the paper is organized as follows: In Section 2
we formulate the HSI classification problem from a Bayesian
perspective. In Section 3 we describe our proposed approach.
In Section 4 we conduct experiments on one synthetic dataset
and two benchmark real HSI datasets and compare with stateof-the-art methods. We conclude in Section 5.

II. P ROBLEM F ORMULATION
Before formulating the HSI classification problem, we first
define the problem-related notations used throughout the paper.
1) Problem notation: Let the HSI dataset be H ∈ Rh×w×d ,
where h and w are the height and width of the spatial
dimensions, respectively, and d is the number of spectral
bands. The set of class labels is defined as K = {1, 2, . . . , K},
where K is the number of classes for the given HSI dataset.
The set of all patches extracted from HSI data H is denoted
as X = {x1 , x2 , . . . , xn }, where xi ∈ Rk×k×d , k is the patch
size in the spatial dimension, n = hw represents the total
number of extracted patches (sliding step is 1 and padding is
used). Each patch xi is first fed into a CNN (shown in Fig.1)
and then a flattened vector zi is output. The corresponding
label of zi is yi ∈ K, which is the label of the spectral
vector corresponding to the center of xi . For the purpose
of simplicity, we denote (xi , yi ) as a sample, which doesn’t
mean all pixel vectors in xi belong to the same class yi
but means that the extracted spatial-spectral feature vector zi
belongs to class yi . In other words, we use the 3-D patches
only to obtain spatial-spectral feature for the corresponding
central pixel vector of this patch. The label set is defined as
y = {y1 , y2 , . . . , yn }. We define the training set for class k
(k)
as Dl(k) = {(x1 , y1 ), . . . , (xl(k) , yl(k) )} and thus the entire
(1)
(2)
(K)
training set can be denoted as Dl = {Dl(1) , Dl(2) , . . . , Dl(K) },
PK (k)
where l =
 n is the total number of training
k=1 l
patches.
2) Problem setup: The goal of HSI classification is to
assign a label yi to the central pixel vector of each patch
xi , i = 1, 2, . . . , n. In the discriminative classification framework, the estimation of labels y for observations X can be
obtained by maximizing a distribution P(y|X , Θ), where Θ
is the parameters of classifier. In our framework, we seek to
combine the spatial modeling power of a Markov random field
with the discriminative power of deep learning. Therefore, we

3

i=1

e = [e
e2T ; . . . ; y
enT ] ∈ Rn×K is an intermediate variwhere y
y1T ; y
ei ∈ RK provides an initial probabilistic
able, in which each y
label pseudo-annotations for each patch xi . In our work, we
define the distribution P(e
yi |xi , Θ) as
(
ei = f (xi ; Θ)
1, y
(2)
P(e
yi |xi , Θ) =
0, otherwise
e is
where f is a classifier. After the pseudo-annotations y
obtained. Then, a second classifier P(y|e
y) takes this collece and outputs the classification
tion of pseudo-annotations y
b = [b
b =
label y
y1 ; yb2 ; . . . ; ybn ], which can be made as y
arg maxy∈Kn P(y|X , Θ).
e is computationally inHowever, since a sum over all y
tensive, we simplify this in two ways for computational
b by maximizing
convenience. First, instead of only learning y
e ) by
P(y|X , Θ), we learn a point estimate of the pairs (b
y, y
e |X , Θ). A goal is then to solve
maximizing P(y, y
)
(
n
X
e ) = arg max log P(y|e
log P(e
yi |xi , Θ) . (3)
(b
y, y
y) +
y,e
y

i=1

Our second simplification takes the form of breaking the optimization problem in Eq. (3) into two subproblems. First, we
e based on the classifier f .
calculate label pseudo-annotations y
e , the current classification
Then, given the annotation values y
b can be obtained by maximizing log P(y|e
results y
y). Thus Eq.
(3) can be solved by the following steps:
= f (xi ; Θ∗ ), i = 1, 2, . . . , n

ei
y
b
y

=

arg max log P(y|e
y).

(4)
(5)

y∈Kn

Here it should be emphasized that classifier f is first trained
on the training data set Dl and Θ∗ is the learned parameters of
the classifier. After obtaining the current classification labels
b on the whole HSI, next we train the classifier again using all
y
b . The two steps
the patches X and their corresponding labels y
are iteratively implemented until some criteria are satisfied. In
this way, we can make full use of the unlabeled patches and
thus more information can be incorporated in our framework.
In this paper, in order to fully utilize both spatial and
spectral information of HSI, we design a new approach based
on this classification framework. Specifically, we firstly adopt a
convolutional neural network (CNN) as the classifier f , which
can not only help to extract spectral-spatial features for HSI,
ei for each xi .
but also can provide label pseudo-annotations y
Then, we place a smoothness prior on labels y to further
enforce spatial consistency. In the following section, we will
introduce these two aspects in detail.

5x5
@100

2x2

3x3
@300

2x2

200

100

K

kxkxd

Fig. 1. The network structure of the CNN used as our discriminative classifier.

III. P ROPOSED A PPROACH
In this section, we first introduce the convolutional neural
network (CNN) classifier for HSIs. Then, we further consider
spatial information by placing a smoothness prior on the labels
y and formulate the classification task as a labeling problem
in an MRF. Finally, we summarize the proposed classification
algorithm.
A. Discriminative CNN-based classifier
Convolutional neural networks (CNN) play an important
role in many computer vision tasks. They consist of combinations of convolutional layers and fully connected layers.
Compared with the standard fully connected feed-forward
neural network (also called multi-layer perceptrons, MLP),
CNNs exploit spatial correlation by enforcing a local connectivity pattern between neurons of adjacent layers to achieve
better performance in many image processing problem, such
as image denoising [63], image super-resolution [21], image
deraining [24] and image classification [37].
In this paper, we adopt the CNN structure shown in Figure 1.
This network contains one input layer, two pairs of convolution
and max pooling layers, two fully connected layers and one
output layer (The tuning of the network structure will be
discussed in Section IV). The detailed parameter settings in
each layer are also shown in Figure 1. For the hyperspectral
image classification task, each sample is a 3D patch of size
k × k × d. Next, we introduce the flow of processing each
sample patch xi at each layer of the CNN.
We note that there is no need to flatten the sample patch xi
into a 1-dimensional vector before it is fed into the input layer.
Therefore, the input size of the first layer is k × k × d. First,
the sample patch is input into the first convolutional layer,
followed by a max pooling operation. The first convolutional
layer filters the k × k × d sample patch using 100 filters of
size 5 × 5 × d. After this convolution, we have 100 feature
maps each n1 × n1 , where n1 = k − 4. We then perform max
pooling on the combined set of features of size n1 × n1 × 100.
The kernel size of the max pooling layer is 2 × 2 and thus
the pooled feature maps have size of n2 × n2 × 100, where
n2 = dn1 /2e. This set of pooled feature maps are then passed
through a second pair of convolutional and max pooling layers.
The convolutional layer contains 300 filters of size 3×3×100,
which filters the pooled feature maps into new feature maps of
size n3 ×n3 ×300, where n3 = n2 −2. Again these new feature

Class Probability

(1)

Fully Connected Layer

e
y

HSI Patch

P(e
yi |xi , Θ),

Fully Connected Layer

n
Y

P(y|e
y)

Fully Connected Layer

e
y

Max Pooling Layer

=

X

P(y|e
y)P(e
y|X , Θ)

Convolutional Layer

=

Max Pooling Layer

e
y

X

Convolutional Layer

let this distribution be of the form
X
e |X , Θ)
P(y, y
P(y|X , Θ) =

4

maps are input into the second max pooling layer with kernel
size 2 × 2 and turned into a second set of pooled feature maps
of size n4 ×n4 ×300, where n4 = dn3 /2e. Finally, the second
pooled feature maps are flattened into a 1-dimensional vector
xpool2 and input to the fully connected layer. The computation
of the next three fully connected layers are as follows:
f (5) (xpool2 )
f

(6)

f

(7)

(xpool2 )
(xpool2 )

=
=
=

σ(W(5) xpool2 + b(5) ),
σ(W
W

(6) (5)

(6)

f

(xpool2 ) + b

(7) (6)

(7)

f

(xpool2 ) + b

,

),

(6)

B. Label smoothness prior and MRF optimization

(7)

e are obtained using the
After the label pseudo-annotations y
current trained CNN. Then, we begin to present our core model
based on the simplifications from Eq. 1 to Eq. 5. Specifically,
b we need to solve
in order to obtain the classification results y
the following optimization problem

(8)

where W(5) , W(6) and W(7) are weight matrices, b(5) ,
b(6) and b(7) are the biases of the nodes, and σ(.) is the
non-linear activation function. In our network, the activation functions in the convolutional layers and fully connected layers are all selected to be the rectified linear
unit (ReLU) function. In order to simplify the notation,
we denote W = {W(1) , W(3) , W(5) , W(6) , W(7) } and
b = {b(1) , b(3) , b(5) , b(6) , b(7) }, where {W(1) , W(3) } and
{b(1) , b(3) } are the weight matrices and biases in the convolutional layers respectively. Thus the aforementioned classifier
parameters Θ are {W, b} in the CNN.
The final vector f (7) ∈ RK is then passed to the
softmax (normalized exponential) function, which gives a
distribution on the label. As defined previously, the label
ei for a given sample xi . Then the
pseudo-annotations is y
pseudo-classification label can be obtained according to yic =
arg maxk∈K yeik . Therefore, given the training set (Dl in the
b )) where we represent
first updating and then the whole (X , y
each scalar label yi ∈ K of a training sample as a one-hot
vector yi ∈ RK , the cross-entropy (CE) loss function can be
computed as
1X
(W,b)
ei
CE(yi , y
),
E(W, b) =
l i
K

1 XX
(W,b)
eik
yik log y
,
= −
l i

(9)

k=1

where W and b is the parameter set defined above and here
(W,b)
ei as y
ei
we denote y
in order to emphasize that the pseudoei is based on the learned CNN. The optimization
annotations y
of the loss function Eq. (9) is conducted by using the stochastic
gradient descent (SGD) algorithm. In the tth iteration, the
weight and bias are updated by
∂E(W, b)
|Wt ,
(10)
∂W
∂E(W, b)
bt+1 = bt − α
|bt ,
(11)
∂b
where the gradients with respect to W and b are calculated
using the back-propagation algorithm [52] and α is the learning rate, which is set as 0.001 in our experiments. Here it
should be noted that unlike some CNN applications, where
the CNN parameters are first pre-trained on some samples
prepared in advance and then fine-tuned on the new dataset,
our proposed CNN-MRF method is directly applied to each
dataset without pre-training. In all our experiments, the CNN
weight parameters W are initialized by random standard
Wt+1

=

Wt − α

normal distribution and the bias parameters b are initialized
with zeros. In order to alleviate the issue of overfitting of the
proposed CNN in the training phase, dropout strategy [56] is
also adopted in the 5th and 6th fully connected layers and the
dropout rate is set as 0.5 in our experiments.

b
y

=

arg max log P(y|e
y)
y∈Kn

=

arg max log P(e
y|y) + log P(y)
y∈Kn

=

arg max
y∈Kn

n
X

log P(e
yi |yi ) + log P(y)

(12)

i=1

e can be regarded as new features for this problem,
where y
log P(e
yi |yi ) is the log-likelihood and P(y) is the label prior.
Specifically, in our work, the log-likelihood log P(e
yi |yi ) is
defined as
log P(e
yi |yi ) =

K
X

1{yi = k} log yeik ,

(13)

k=1

where 1{·} is an indicator function.
In image segmentation tasks, it is probable that adjacent
pixels have the same label. The exploitation of this naive prior
information can often dramatically improve the segmentation
performance. In this paper, we enforce a spatially smooth prior
on labels y to encourage neighboring pixels to belong to the
same class. This smoothness prior distribution on labels y is
defined as
1 Pn P
(14)
P(y) = eµ i=1 j∈N (i) δ(yi −yj ) ,
Z
where Z is a normalization constant for the distribution, µ
is the label smoothness parameter, N (i) is the neighboring
pixels of pixel i and δ(·) is a function defined as: δ(0) = 1
and δ(y) = −1 for y 6= 0. We note that the pairwise
interaction terms δ(yi − yj ) obtain higher probability when
neighboring labels are equal than when they are not equal.
In this way, this smoothness prior can encourage piecewise
smooth segmentations.
Based on Eq. (13) and Eq. (14), the final classification
model (12) is thus given by


n X
K
n X
X

X
b = arg max
y
1{yi = k}log yeik +µ
δ(yi −yj ) .

y∈Kn 
i=1 k=1

i=1j∈N (i)

(15)
This objective function contains many pairwise interaction
terms and is a challenging combinatorial optimization problem. It can also be regarded as an MRF model in which an
undirected model graph G =< V, E > is defined on the whole
image, where graph node sets V correspond to pixels, and the
undirected edge set E represents the neighboring relationship

5

between the pixels [9]. We define a random variable yi on each
node vi ∈ V, thus allowing the labels y to form a Markov
random field. The objective function is its energy function,
of which the first term represents the cost of a pixel being
assigned with different classes. The larger the probability of a
pixel belonging to a certain class, the more probable that the
pixel is assigned the corresponding label. The second term
of the energy function encourages the labels of neighboring
pixels to be the same.
This optimization problem of the MRF is NP-hard. Many
approximating algorithms have been proposed, such as the
graph cut method [8], [36], belief propagation [65] and message passing [35]. In this paper, we use the α-expansion mincut method [8] because of its good performance and fast
computation speed.

Training
data 𝓓𝓓𝒍𝒍
Train CNN

� using
Update 𝒚𝒚
Graph Cut

Algorithm 1 CNN-MRF classification algorithm for HSI
Input: HSI patches X , training data Dl , learning rate α,
smoothness parameter µ, batch size nbatch .
b.
Output: Labels y
1: Train the CNN classifier using Dl ;
e on X and update label y
b = α-Expansion(e
2: Compute y
y,µ).
b };
3: Train the CNN classifier using {X , y
e on X and update label y
b = α-Expansion(e
4: Compute y
y,µ).
5: Repeat step 3 and 4 until the stopping criterion is satisfied.

IV. E XPERIMENTS
To test the effectiveness of our proposed CNN-MRF classification algorithm in different scenarios, we conduct experiments on one synthetic dataset and two real-world benchmark
datasets.1 For comparison, we consider several state-of-the-art
HSI classification methods, including support vector machine
graph-cut method (SVM-GC), subspace multinomial logistic
regression with multilevel logistic prior (MLRsubMLL)2 [38]
and support vector machine based on the 3-dimensional discrete wavelet transform method (SVM-3DDWT-GC)3 [12].
It should be noted that the competing methods SVM-GC
and MLRsubMLL are methods that integrate the MRF into
the SVM and MLRsub methods, respectively. In order to
distinguish the classification methods with MRF and without
MRF, we denote them as classification methods and regularized classification methods, respectively. Besides, we also
compare with three deep learning methods: SS-DCNN [67],
SPP-DCNN [66] and DC-CNN [68]. Our proposed CNN-MRF
approach is implemented in Python using the Tensorflow [1]
1 Code

is available at: https://github.com/xiangyongcao/CNN HSIC MRF.
is available at: https://www.lx.it.pt/∼jun/demos.html
3 Code is available at: https://github.com/xiangyongcao/3DDWT-SVM-GC
2 Code

Data 𝓧𝓧

Train CNN

� using
Update 𝒚𝒚
Graph Cut

C. CNN-MRF classification algorithm
After introducing how to update CNN’s parameters and
compute class labels respectively, we can summarize our
CNN-MRF classification algorithm in detail in Algorithm 1.
Besides, the corresponding flowchart of this algorithm is also
depicted in Figure 2.

Data 𝓧𝓧

�
Output 𝒚𝒚

Fig. 2. The flowchart of the proposed CNN-MRF algorithm.

library on a server with Nvidia GeForce GTX 1080 and Tesla
K40c. All non-deep algorithms are run in Matlab R2014b. For
comparisons with other deep learning models, we use the best
reported results for these algorithms.
All methods are compared numerically using the following
three criteria [12]: overall accuracy (OA), average accuracy
(AA) and the kappa coefficient (κ). OA represents the number
of correctly classified samples divided by the total number
of test samples, AA denotes the average of individual class
accuracies, and κ involves both omission and commission
errors and gives a good representation of the the overall
performance of the classifier. For all the three criteria, a larger
value indicates a better classification performance.
A. Synthetic HSI data
In this section, we first generate a synthetic HSI dataset.
Then, we use this dataset to evaluate in detail the sensitivity
of performance to different parameters setting of our CNN
structure. Finally, we compare our method with other competing methods on this dataset using the tunning CNN structure.
1) Generation of Synthetic HSI data: To generate a
synthetic HSI data, five endmembers are first extracted
randomly from a real scene with 162 bands in ranges
400−2500 nm, and then 40000 vectors are generated as a
sum of Gaussian fields with constraints so as to respect
the abundance-nonnegative-constraint (ANC) and abundancesum-to-one-constraint (ASC). Finally, this dataset is generated
using a Generalized Bilinear Mixing Model (GBM) [69]:
z =

K
X
i=1

ai ei +

K−1
X

K
X

γij ai aj ei

ej + n, (16)

i=1 j=i+1

with class number K = 5. Here z is the simulated pixel vector,
γij are selected uniformly from [0, 1], ei , i = 1, . . . , 5 are the

6

five end-members,
n is the Gaussian noise with an SNR of 30
PK
dB, ai ≥ 0 and i=1 ai = 1.
2) Impact of parameter settings: In this section, we evaluate in detail the sensitivity of performance to different
parameters settings of our CNN structure. In the following
experiments, we use the synthetic dataset and randomly choose
1% training samples from each class as training data and the
remaining for testing.
(1) Kernel Size: First, we test the impact of different kernel
sizes. We fix the kernel size of the second layer as 3 and
evaluate the performance by varying the kernel size in the
first layer. Table I shows these results. From this table we can
conclude that larger kernel sizes can obtain better results. This
is because more structure and texture can be captured using a
larger kernel size. Therefore, in our experiments, we set the
kernel size of the first layer as 5.
TABLE I
OVERALL ACCURACY (%) WITH DIFFERENT KERNEL SIZES .

Kernel Size
OA

1
96.07

2
96.71

3
97.10

4
98.26

5
99.51

(2) Network Width: Next, we evaluate the impact of network
width in the convolutional layer on the classification results.
Fixing the network width of the first convolutional layer to
100, we test the performance by changing the width of the
second convolutional layer. These results are displayed in
Table II. It can be observed that the results are not very
sensitive to the network width and thus in our experiments,
we select 200 as the default setting of the network width in
the second convolutional layer.
TABLE II
OVERALL ACCURACY (%) WITH DIFFERENT NETWORK WIDTHS .

Network Width
OA

50
99.27

100
99.36

200
99.51

300
99.32

500
99.35

(3) Network Depth: We also conduct experiments to test
the impact of different network depths by reducing or adding
non-linear layers. We train and test on 5 networks with depths
5, 7, 9, 11 and 13. The experimental results are summarized
in Table III. From this table it can be seen that increasing
the network depth does not always generate better results,
which may be a result of the gradient vanishing problem [28]
encountered by deeper neural networks. For this experiment,
the best performance can be achieved by setting the network
depth to 7, which we use as the default setting for the network
depth in our experiments.
TABLE III
OVERALL ACCURACY (%) WITH DIFFERENT NETWORK DEPTHS .

Network Depth
OA

5
99.35

7
99.52

9
99.48

11
99.45

13
99.41

(4) Patch Size: We also investigate the performances of our
proposed method with respect to different data patch sizes k =

{1, 3, 5, 9, 13}. These results are shown in Table IV. As can be
seen, a larger patch size of the sample generates better results.
This is because more spatial information is incorporated into
the training process. However, it also takes more computation
time to train the network with an increasing patch size. Thus
we choose 9 as the default setting of patch size as a tradeoff
between performance and the running time.
TABLE IV
OVERALL ACCURACY (%) WITH DIFFERENT PATCH SIZES .

Patch Size
OA

1
95.12

3
96.40

5
98.13

9
99.43

13
99.52

(5) Smoothness parameter: Finally, we test the impact of
different smoothness parameter µ on the obtained classification
results. Table V shows these results, where it can be observed
that setting µ = 20 produces the best classification results.
Thus we use this value in our experiments.
TABLE V
OVERALL ACCURACY (%) WITH DIFFERENT SMOOTH PARAMETERS .

Smooth Parameter
OA

1
98.43

3
98.61

5
98.70

10
98.89

20
99.53

30
99.31

50
99.35

For other parameters, such as the kernel size in max pooling
layer which is set as 2, the number of nodes in the two fully
connected layers, which are set as 200 and 100 respectively, we
observe consistent performance under variations. Therefore,
we didn’t report their impact to the performance in this
section. The learning rate α is set as 0.001 and the batch
size nbatch is set as 100 for synthetic and Indian Pines dataset
since consistent performance is also observed under variations.
While for Pavia University dataset the learning rate is set
to 0.0001 and the batch size nbatch is set as 150. For the
following experiments, we adopt the same network structure.
3) Experimental result of the synthetic HSI data: Using the
parameter settings above for the CNN structure, we compare
our CNN-MRF method with other competing methods. The
final classification results are shown in Table VI. From Table
VI, we see that the proposed CNN-MRF method achieves
better performance with respect to OA, AA and κ than other
methods. We emphasize that all non-deep learning methods
(namely SVM-GC, MLRsubMLL and SVM-3DDWT-GC) in
this experiment use the MRF to post-process the classification
map, and the differences between them and our method are
the way to extract features and the iteratively updating the
CNN parameters and class labels. Thus the improvement in our
method here is due to the use of CNN and the integration of
CNN with MRF. Comparing with other deep-learning method,
our CNN-MRF method also outperforms them. For better
visualization, we also demonstrate the final classification map
in Figure 3. From Figure 3, we can see that the classification
maps obtained by our method are visually closer to the ground
truth map than the other methods, which is consistent with the
quantitative metric OA.

7

TABLE VI
OVERALL ACCURACIES (%), AVERAGE ACCURACIES (%), KAPPA STATISTICS AND RUNNING TIME OF ALL COMPETING METHODS ON THE SYNTHETIC
DATASET.
Class
OA
AA
κ

SVM-GC [60]
96.44
94.55
95.23

MLRsubMLL [38]
97.92
95.74
96.03

SVM-3DDWT-GC [12]
93.89
91.37
92.84

SS-DCNN [67]
94.75
92.46
93.57

SPP-DCNN [66]
95.63
93.24
94.46

DC-CNN [68]
98.78
96.59
97.21

CNN-MRF
99.55
97.83
98.26

Ground-truth

SVM-GC (96.44%)

MLRsubMLL(97.92%)

SVM-3DDWT-GC(93.89%)

SS-DCNN (94.75%)

SPP-DCNN (95.63%)

DC-CNN (98.78%)

CNN-MRF(99.55%)

Fig. 3. Classification maps obtained by all competing methods on the synthetic dataset (overall accuracies are reported in parentheses).

TABLE VII
S TATISTICS OF THE I NDIAN P INES DATA SET, INCLUDING THE NAME , THE
NUMBER OF TRAINING , TEST AND TOTAL SAMPLES FOR EACH CLASS .

No
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

Class
Name
Alfalfa
Corn-no till
Corn-min till
Corn
Grass-pasture
Grass-trees
Grass-pasture-mowed
Hay-windrowed
Oat
Soybean-no till
Soybean-min till
Soybean-clean
Wheat
Woods
Buildings-Grass-Trees-Drives
Stone-Steel-Towers
Total

Train
5
143
83
24
49
73
3
48
2
98
246
60
21
127
39
10
1025

Samples
Test
41
1285
747
213
434
657
25
430
18
847
2209
533
184
1138
347
83
9224

Total
46
1428
830
237
483
730
28
478
20
972
2455
593
205
1265
386
93
10249

B. Real HSI Data
We next test our method on two real HSI benchmark
datasets, the Indian Pines dataset and Pavia University dataset,
in a variety of experimental settings.

1) AVIRIS Indian Pines Data: This data set was gathered by
the Airborne Visible/Infrared Imaging Spectrometer (AVIRIS)
sensor over the Indian Pines test site in North-western Indiana
in June 1992. The original dataset contains 220 spectral reflectance bands in the wavelength range 0.4−2.5 µm, of which
20 bands cover the region of water absorption. Unlike other
methods which remove the 20 polluted bands, we keep all 220
bands in our experiments. This HSI has a spectral resolution
of 10 nm and a spatial resolution of 20 m by pixel, and the
spatial dimension is 145×145. The ground truth contains 16
land cover classes. This dataset poses a challenging problem
because of the significant presence of mixed pixels in all
available classes and also because of the unbalanced number
of available labeled pixels per class. A sample band of this
dataset and the related ground truth categorization map are
shown in Figure 4.
To first evaluate our proposed CNN-MRF method in the
scenario of limited training samples, we randomly choose
10% of the available labeled samples for each class from
the reference data, which is an imbalanced training sample
case, and use the remaining samples in each class for testing.
The training and testing sets are summarized in Table VII.
This experiment is repeated 20 times for each method and the
average performance is reported.
Additionally, some parameters need to be set in advance in
these experiments. For the competing methods, their parame-

8

TABLE VIII
I NDIVIDUAL CLASS , OVERALL , AVERAGE ACCURACIES (%) AND KAPPA STATISTICS OF ALL METHODS ON THE I NDIAN P INES IMAGE TEST SET.

Class
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
OA
AA
Kappa

SVM
73.17
62.65
52.88
32.39
91.24
92.09
36.00
95.58
0
61.44
86.92
76.36
91.85
97.01
48.13
91.57
77.02
68.08
73.49

(a)

Classification algorithms
MLRsub SVM-3DDWT
46.34
63.41
40.93
89.81
26.24
91.97
17.37
80.75
70.97
96.77
94.37
98.78
18.18
56.00
96.51
100
22.22
94.44
25.06
88.44
78.23
95.79
16.51
94.93
93.48
94.57
99.38
97.72
4.32
79.83
77.11
75.90
63.12
93.19
50.57
87.44
53.13
92.22

CNN
85.23
90.17
93.43
84.19
98.76
99.43
74.52
98.74
100
91.71
95.59
89.97
98.64
97.88
89.95
96.83
94.05
92.97
92.86

SVM-GC [60]
95.12
68.48
56.49
77.00
94.47
97.72
34.42
100
0
75.06
95.47
99.44
98.37
97.45
76.66
98.80
85.92
76.91
83.78

(b)

Fig. 4. Indian Pines image and related ground truth categorization information. (a) The original HSI. (b) The ground truth categorization map (This
figure is better seen by zooming on a computer screen.)

ters are set as their papers suggest: For SVM-based methods,
the RBF kernel parameter γ and the penalty parameter C are
tuned through 5-fold cross validation (γ = 2−8 , 2−7 , . . . , 28 ,
C = 2−8 , 2−7 , . . . , 28 ). For SVM-GC and SVM-3DDWT-GC,
the spatial smoothness parameter β is set as 0.75 advised
by [12], [61]. For MLRsubMLL, the smoothness parameter
µ and the threshold parameter τ are set following [38]. For
our proposed CNN-MRF method, we adopt the same network
structure as the synthetic HSI data.
For the above experimental settings, in order to measure the
performance improvement due to including spatial contextual
information with the MRF, we also report the classification
results of each method without using the MRF prior, and
call these methods SVM, MLRsub, SVM-3DDWT and CNN,
respectively. Classification maps on the Indian dataset are
shown for all methods in Figure 5, and the accuracies (i.e.
individual class accuracy, OA, AA and κ) are reported in Table
VIII.
From Table VIII, we highlight two main results. First,
CNN-MRF and CNN, which are both deep learning methods,
achieve the first and third best performance in terms of the

Regularized Classification algorithms
MLRsubMLL [38] SVM-3DDWT-GC [12]
95.12
82.93
50.04
95.56
13.12
95.72
15.02
95.77
73.04
96.54
98.93
99.39
37.25
0
100
100
0
88.89
19.68
92.22
88.82
98.05
16.51
98.31
99.46
98.37
99.91
99.03
60.52
89.91
83.13
72.29
70.45
94.28
56.82
87.69
65.43
94.85

CNN-MRF
86.52
91.46
96.35
96.22
99.48
99.82
78.00
98.84
100
94.26
96.48
91.85
98.85
98.36
91.54
97.85
96.12
94.75
95.78

three criteria (OA, AA and κ); CNN-MRF attains about 2%
improvement in term of OA in this scenario, with SVM3DDWT-GC trailing slightly behind (94.28%). And CNN falls
behind SVM-3DDWT-GC only 0.23%. This means that using
a CNN for classification without the MRF performs better
than a MRF-based model with other classifier except the
SVM classifier with 3DDWT features. The CNN therefore
significantly helps for this problem. Moreover, as depicted in
Figure 5, the classification maps of the CNN-based methods
are noticeably closer to the ground truth map. Finally, directly
comparing the MRF and non-MRF based methods, we can
conclude that using an MRF prior significantly improves the
classification accuracy of any particular classifier because it
further embeds the spatial smoothness information into the
segmentation stage. Therefore, the superior performance of our
proposed CNN-MRF method can be explained by using the
CNN and MRF strategies simultaneously to fully exploit the
spectral and spatial information in a HSI.
2) ROSIS Pavia University Data: We perform similar experiments on a second real HSI dataset. This HSI was acquired by the Reflective Optics System Imaging Spectrometer
(ROSIS) over the urban area of the University of Pavia in
northern Italy on July 8, 2002. The original dataset consists
of 115 spectral bands ranging from 0.43 to 0.86 µm, of which
12 noisy bands are removed and only 103 bands are retained
in our experiments. The scene has a spatial resolution of 1.3
m per pixel, and the spatial dimension is 610×340. There are
9 land cover classes in this scene and the number of each
class is displayed in Table IX. We show a sample band and
the corresponding ground truth class map in Figure 6.
To evaluate the performance of our proposed CNN-MRF
method using only a small number of labeled training samples,
we randomly chose 40 samples for each class from the ground
truth data for training, which gives a balanced training sample,
and the remaining samples are used for testing. The related

9

SVM (77.02%)

MLRsub(63.12%)

SVM-3DDWT(93.19%)

CNN(94.05%)

SVM-GC (85.92%)

MLRsubMLL(70.45%)

SVM-3DDWT-GC(94.28%)

CNN-MRF(96.12%)

Ground-truth

Fig. 5. Classification maps obtained by all methods on the Indian Pines dataset (overall accuracies are reported in parentheses).

TABLE IX
S TATISTICS OF THE PAVIA U NIVERSITY DATA SET, INCLUDING THE NAME ,
THE NUMBER OF TRAINING , TEST AND TOTAL SAMPLES FOR EACH CLASS .

No
1
2
3
4
5
6
7
8
9

Class
Name
Asphalt
Meadows
Gravel
Trees
Painted metal sheets
Bare Soil
Bitumen
Self-Blocking Bricks
Shadows
Total

Train
40
40
40
40
40
40
40
40
40
360

Samples
Test
6591
18609
2059
3024
1305
4989
1290
3642
907
42416

Total
6631
18649
2099
3064
1345
5029
1330
3682
947
42776

statistics are also summarized in Table IX. As previously,
we repeat this experiment 20 times for each method and
report the average performance. All parameters involved in the
compared methods are tuned in the same way as the previous
Indian Pines experiment. The network structure settings of the
CNN are also the same. Classification and segmentation maps
obtained by all methods on this dataset are illustrated in Figure
7 and accuracies (i.e. individual class accuracy, OA, AA and
κ) are summarized in Table X.
From Table X, we can conclude that, for this dataset, our
CNN-MRF approach again achieves the best performance in
terms of the three quantitative criteria. It is also worth noting
that the CNN without MRF again performs third best with
respect to OA and κ, which means that the CNN plays an important role in improving the classification accuracy. By using
the MRF prior, CNN-MRF obtains about 2% improvement in
terms of the OA compared with the CNN. However, for other
classification algorithms we observe that the MRF can greatly
boost classification accuracy. For example, MLRsubMLL has
about 30% improvement according to OA compared with
MLRsub. SVM-3DDWT-GC has the second best classification OA (94.12%) due to the 3D discrete wavelet transform
(3DDWT), as has been previously studied in [12]. Meanwhile,
it can be seen from Figure 7 that CNN-MRF obtains much
smoother classification map than other methods, which is
consistent with the results shown in Table X. Consequently,
improvement of our CNN-MRF approach can be explained by
the use of CNN and MRF models simultaneously.
C. Limited training data scenarios

(a)

(b)

Fig. 6. Pavia University image and related ground truth categorization
information. (a) The original HSI. (b) The ground truth categorization map.
(This figure is better seen by zooming on a computer screen.)

To analyze the sensitivity of the proposed method to training
sets consisting of limited training samples, we conduct additional experiments in which 0.1%, 0.2%, 0.3%, 0.4% and 0.5%
of each class are randomly selected from the Pavia University
data as training samples and the remaining are used for
testing. For this experiment, we adopt the data augmentation
technique [37] to help the training process of the CNN. The

10

TABLE X
I NDIVIDUAL CLASS , OVERALL , AVERAGE ACCURACIES (%) AND KAPPA STATISTICS OF ALL COMPETING METHODS ON THE PAVIA U NIVERSITY IMAGE
TEST SET.

Class
1
2
3
4
5
6
7
8
9
OA
AA
Kappa

SVM
72.17
64.85
75.57
89.62
97.78
72.42
86.82
67.13
97.57
73.41
79.36
66.23

Classification algorithms
MLRsub SVM-3DDWT
62.10
84.78
46.57
93.62
58.38
85.04
89.75
95.90
99.54
98.77
73.80
94.29
80.78
96.59
66.56
79.82
99.56
100
61.36
91.27
75.23
92.09
53.26
88.55

CNN
96.82
96.80
86.16
98.54
99.88
90.40
86.92
90.94
97.53
94.82
93.77
93.89

SVM (73.41%)

SVM-GC [60]
97.74
67.92
90.38
90.34
99.85
94.75
71.16
68.51
99.67
80.21
86.70
75.09

Regularized classification algorithms
MLRsubMLL [38] SVM-3DDWT-GC [12]
87.07
92.03
96.97
97.33
77.27
86.89
83.90
97.09
99.54
98.70
99.40
98.12
94.50
98.84
64.83
84.46
99.78
100
91.13
94.12
89.25
94.83
88.19
93.55

MLRsub(61.36%)

CNN-MRF
98.02
97.78
88.47
99.17
99.90
93.00
87.47
91.66
98.03
96.18
94.83
94.62

SVM-3DDWT(91.27%)

CNN(94.82%)

SVM-3DDWT-GC(94.12%)

CNN-MRF(96.18%)

Ground-truth

SVM-GC (80.21%)

MLRsubMLL(91.13%)

Fig. 7. Classification maps obtained by all competing methods on the Pavia University dataset (overall accuracies are reported in parentheses).

OA of all methods are displayed in Figure 8. From Figure
8(a), we observe that the classification results of the CNN
method outperform the other methods for each training set
size. Additionally, when the spatial prior is considered using
the MRF, the classification results in Figure 8(b) significantly
improve the corresponding classification results in Figure 8(a),
again indicating that the MRF is an important factor for
improving the classification accuracy.
D. Study on the interaction between CNN and MRF
Furthermore, in order to clearly illustrate the interactive
influence between CNN and MRF, we also report the OA of

our method on the three datasets as a function of iteration. The
results are shown in Table XI. Specifically, in our experiments,
we first train our CNN for 30 epochs using the training data
b every 10 epochs. From
Dl . Then, we update the class labels y
Table XI, it can be seen that the OA first increases quickly
in the first 30 epochs. Then, the OA increases slowly for
about 40 epochs. Finally, the OA decreases a little or has a
slight fluctuation. Therefore, we conclude that the interaction
between the CNN and MRF can help the final classification
OA. In our experiments, we set the maximum epoch as 60 and
report the output as the results of our method.

95

95

90

90
Overall accuracy (%)

Overall accuracy (%)

11

85
80
75

SVM
MLRsub
SVM-3DDWT

70

85
80
75

SVM-GC
MLRsubMLL
SVM-3DDWT-GC
CNN-MRF

70

CNN

65
0.1

0.2
0.3
0.4
Percentage of training samples per class(%)

0.5

65
0.1

0.2
0.3
0.4
Percentage of training samples per class(%)

0.5

Fig. 8. Overall accuracy (%) obtained by all competing methods with different proportions of training samples on Pavia University dataset. (a) Classification
results. (b) Regularized classification results.

TABLE XI
T HE CHANGE OF OVERALL ACCURACY (%) AS EPOCH GOES .
Epoch
Synthetic
Indian Pines
PaviaU

10
92.53
87.29
83.71

20
96.72
92.67
92.71

30
98.85
94.05
94.82

40
99.44
95.36
95.71

50
99.52
95.91
95.92

60
99.58
96.12
96.18

70
99.44
95.93
96.17

80
99.30
96.03
96.18

90
99.12
95.98
96.16

TABLE XII
OVERALL ACCURACY, AVERAGE ACCURACY AND KAPPA COEFFICIENT (%)
OF ALL THE DEEP LEARNING METHODS ON THE I NDIAN P INES DATASET.
100
99.25
95.97
96.18

E. Comparison with other deep learning methods
In order to further evaluate the performance of our proposed CNN-MRF method, we compare with three recent deep
learning HSI classification algorithms: SS-DCNN [67], SPPDCNN [66] and DC-CNN [68]. We use the two previous real
datasets for comparison. For fair comparison, we choose the
training samples with same proportion as was done in [68].
Specifically, for the Indian Pines dataset, we choose 10%
samples from each class as training set, and for the Pavia
University dataset we choose 5% samples from each class as
training data. We also apply the data augmentation strategy
in the two experiments since all deep learning algorithms can
adopt this strategy. The results for each algorithm are displayed
in Table XII and Table XIII.4 Results for Indian Pines are
shown in Table XII, where it can be seen that our proposed
CNN-MRF method achieves improved performance compared
with other deep methods in terms of OA, AA and κ. We
show results for the Pavia University data in Table XIII. We
again observe that our method achieves the best performance
compared with other methods. From Table XII and XIII, we
can conclude that the running time of our proposed method
is a potential drawback. Here the running time reported in
this experiment includes both training and testing time. (This
could be further addressed with high-performance computing
resources and techniques.)
F. Comparison with other label regularization methods
In order to evaluate the performance of MRF, we compare
it with other label regularization methods, such as median

OA
AA
kappa(κ)
Time

SS-DCNN [67]
90.76
85.52
89.44
255.32

SPP-DCNN [66]
91.60
93.96
90.43
328.18

DC-CNN [68]
98.76
98.50
98.58
4860.57

CNN-MRF
99.32
99.27
99.21
1454.62

TABLE XIII
OVERALL ACCURACY, AVERAGE ACCURACY AND KAPPA COEFFICIENT (%)
OF ALL THE DEEP LEARNING METHODS ON THE PAVIA U NIVERSITY
DATASET.

OA
AA
kappa(κ)
Time

SS-DCNN [67]
93.34
92.20
91.95
266.11

SPP-DCNN [66]
94.88
93.29
93.21
333.90

DC-CNN [68]
99.68
99.50
99.58
3421.25

CNN-MRF
99.71
99.55
99.63
1156.43

filter [30] and majority voting [7]. The other two methods
are implemented just by replacing the MRF optimization
with median filter method and majority voting method in our
algorithm framework (we call the two methods as CNN-MF
and CNN-MV respectively). All the methods are compared
on the above synthetic data, Indian Pines data and Pavia
University data with the same experimental settings as shown
earlier. For median filter and majority voting methods, we set
the window size as {3,5,7} and report the best result. The
experimental results are shown in Table XIV and Figure 9.
From Table XIV and Figure 9, we can easily observe that
in MRF obtains the best performance comparing with other
two simple label regularization methods on all the datasets.
Thus we adopt MRF as the label regularization method in our
algorithm.
V. C ONCLUSIONS

4 The

results for SS-DCNN, SPP-DCNN and DC-CNN are taken from [68].
The running time of Table VIII and Table XII is different since data
augmentation strategy is adopted in this experiment of Table XII.

In this paper, we proposed a novel technique for HSI
classification that incorporates both spectral and spatial infor-

12

Ground-truth

CNN(98.89%)

CNN-MF(93.92%)

CNN-MV(94.10%)

CNN-MRF(99.25%)

Ground-truth

CNN(94.05%)

CNN-MF(95.19%)

CNN-MV(95.37%)

CNN-MRF(95.52%)

Ground-truth

CNN(94.82%)

CNN-MF(95.44%)

CNN-MV(95.55%)

CNN-MRF(95.69%)

Fig. 9. Classification maps obtained by all competing methods on the Pavia University dataset (overall accuracies are reported in parentheses).

TABLE XIV
OVERALL ACCURACY (%) OF ALL THE LABEL REGULARIZATION
METHODS ON THE SYNTHETIC , I NDIAN PINES AND PAVIA U NIVERSITY
DATASETS .

Synthetic
Indian Pines
PaviaU

CNN
98.89
94.05
94.82

CNN+Median Filter
93.92
95.19
95.44

CNN+Majority Voting
94.10
95.37
95.55

CNN-MRF
99.25
95.52
95.69

mation in a unified Bayesian framework. Specifically, we use
a convolutional neural network (CNN) in combination with
a Markov random field to classify HSI pixel vectors in a
way fully takes spatial and spectral information into account.
We then efficiently learn the classfication result by iteratively
updating the CNN parameters and the class labels of all
HSI pixel vectors. Experimental results on one synthetic HSI
dataset and two real benchmark HSI datasets show that our
method outperforms state-of-the-art methods, including deep
and non-deep models. In the future, we will further consider
the HSI classification task in unsupervised settings. Besides,
we will also try to extend our model to the popular deep
generative models, such as variational autoencoder (VAE) [33]
and generative adversarial network (GAN) [26]. Furthermore,

we hope to design more powerful regularization regimes,
extending the employed MRF one, for different application
scenarios of this problem in our future research. We will also
consider dimensionality reduction techniques in our methods
to further make the method more efficiently computed, especially on large-scaled datasets.
R EFERENCES
[1] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S.
Corrado, A. Davis, J. Dean, M. Devin, et al. Tensorflow: Large-scale
machine learning on heterogeneous distributed systems. arXiv preprint
arXiv:1603.04467, 2016.
[2] T. V. Bandos, L. Bruzzone, and G. Camps-Valls. Classification of
hyperspectral images with regularized linear discriminant analysis. IEEE
Transactions on Geoscience and Remote Sensing, 47(3):862–873, 2009.
[3] J. A. Benediktsson, J. A. Palmason, and J. R. Sveinsson. Classification
of hyperspectral data from urban areas based on extended morphological
profiles. IEEE Transactions on Geoscience and Remote Sensing,
43(3):480–491, 2005.
[4] J. Besag. Spatial interaction and the statistical analysis of lattice systems.
Journal of the Royal Statistical Society. Series B (Methodological), pages
192–236, 1974.
[5] M. Bhatt and U. B. Desai. Robust image restoration algorithm using
markov random field model. Graphical Models and Image Processing
(CVGIP), 56(1):61–74, 1994.
[6] H. Bischof and A. Leonardis. Finding optimal neural networks for
land use classification. IEEE transactions on Geoscience and Remote
Sensing, 36(1):337–341, 1998.
[7] R. S. Boyer and J. S. Moore. MJRTYA Fast Majority Vote Algorithm.
Springer Netherlands, 1991.

13

[8] Y. Boykov, O. Veksler, and R. Zabih. Fast approximate energy minimization via graph cuts. IEEE Transactions on Pattern Analysis and
Machine Intelligence, 23(11):1222–1239, 2001.
[9] Y. Y. Boykov and M.-P. Jolly. Interactive graph cuts for optimal
boundary & region segmentation of objects in nd images. In IEEE
International Conference on Computer Vision (ICCV), 2001, volume 1,
pages 105–112. IEEE, 2001.
[10] G. Camps-Valls and L. Bruzzone. Kernel-based methods for hyperspectral image classification. IEEE Transactions on Geoscience and Remote
Sensing, 43(6):1351–1362, 2005.
[11] G. Camps-Valls, L. Gómez-Chova, J. Calpe-Maravilla, J. D. Martı́nGuerrero, E. Soria-Olivas, L. Alonso-Chordá, and J. Moreno. Robust
support vector method for hyperspectral data classification and knowledge discovery. IEEE Transactions on Geoscience and Remote sensing,
42(7):1530–1542, 2004.
[12] X. Cao, L. Xu, D. Meng, Q. Zhao, and Z. Xu. Integration of 3dimensional discrete wavelet transform and markov random field for
hyperspectral image classification. Neurocomputing, 226:90–100, 2017.
[13] L.-C. Chen, G. Papandreou, I. Kokkinos, K. Murphy, and A. L. Yuille.
Semantic image segmentation with deep convolutional nets and fully
connected crfs. arXiv preprint arXiv:1412.7062, 2014.
[14] Y. Chen, Z. Lin, X. Zhao, G. Wang, and Y. Gu. Deep learning-based
classification of hyperspectral data. IEEE Journal of Selected topics in
applied earth observations and remote sensing, 7(6):2094–2107, 2014.
[15] Y. Chen, N. M. Nasrabadi, and T. D. Tran. Hyperspectral image classification using dictionary-based sparse representation. IEEE Transactions
on Geoscience and Remote Sensing, 49(10):3973–3985, 2011.
[16] Y. Chen, N. M. Nasrabadi, and T. D. Tran. Hyperspectral image
classification via kernel sparse representation. IEEE Transactions on
Geoscience and Remote Sensing, 51(1):217–231, 2013.
[17] Y. Chen, X. Zhao, and X. Jia. Spectral–spatial classification of hyperspectral data based on deep belief network. IEEE Journal of Selected
Topics in Applied Earth Observations and Remote Sensing, 8(6):2381–
2392, 2015.
[18] D. L. Civco. Artificial neural networks for land-cover classification and
mapping. International Journal of Geographical Information Science,
7(2):173–186, 1993.
[19] M. Crawford and W. Kim. Manifold learning for multi-classifier systems
via ensembles. Multiple Classifier Systems, pages 519–528, 2009.
[20] M. Dalla Mura, A. Villa, J. A. Benediktsson, J. Chanussot, and
L. Bruzzone. Classification of hyperspectral images by using extended
morphological attribute profiles and independent component analysis.
IEEE Geoscience and Remote Sensing Letters, 8(3):542–546, 2011.
[21] C. Dong, C. C. Loy, K. He, and X. Tang. Image super-resolution using
deep convolutional networks. IEEE transactions on pattern analysis and
machine intelligence, 38(2):295–307, 2016.
[22] J. M. Duarte-Carvajalino, G. Sapiro, M. Vélez-Reyes, and P. E. Castillo.
Multiscale representation and segmentation of hyperspectral imagery
using geometric partial differential equations and algebraic multigrid
methods. IEEE Transactions on Geoscience and Remote Sensing,
46(8):2418–2434, 2008.
[23] M. Fauvel, Y. Tarabalka, J. A. Benediktsson, J. Chanussot, and J. C.
Tilton. Advances in spectral-spatial classification of hyperspectral
images. Proceedings of the IEEE, 101(3):652–675, 2013.
[24] X. Fu, J. Huang, X. Ding, Y. Liao, and J. Paisley. Clearing the skies:
A deep network architecture for single-image rain removal. IEEE
Transactions on Image Processing, 26(6):2944–2956, 2017.
[25] I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press,
2016. http://www.deeplearningbook.org.
[26] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley,
S. Ozair, A. Courville, and Y. Bengio. Generative adversarial nets. In
Advances in neural information processing systems, pages 2672–2680,
2014.
[27] K. He, X. Zhang, S. Ren, and J. Sun. Spatial pyramid pooling in deep
convolutional networks for visual recognition. In European Conference
on Computer Vision (ECCV), 2014, pages 346–361. Springer, 2014.
[28] S. Hochreiter. The vanishing gradient problem during learning recurrent
neural nets and problem solutions. International Journal of Uncertainty,
Fuzziness and Knowledge-Based Systems, 6(02):107–116, 1998.
[29] W. Hu, Y. Huang, L. Wei, F. Zhang, and H. Li. Deep convolutional neural
networks for hyperspectral image classification. Journal of Sensors,
2015, 2015.
[30] T. Huang, G. Yang, and G. Tang. A fast two-dimensional median
filtering algorithm. IEEE Transactions on Acoustics Speech and Signal
Processing, 27(1):13–18, 1979.
[31] S. Ji, W. Xu, M. Yang, and K. Yu. 3d convolutional neural networks for
human action recognition. IEEE transactions on pattern analysis and
machine intelligence, 35(1):221–231, 2013.
[32] S. Jia, X. Zhang, and Q. Li. Spectral–spatial hyperspectral im-

[33]
[34]
[35]
[36]
[37]
[38]

[39]

[40]

[41]

[42]

[43]

[44]
[45]
[46]
[47]

[48]

[49]

[50]

[51]

[52]
[53]

[54]

age classification using regularized low-rank representation and sparse
representation-based graph cuts. IEEE Journal of Selected Topics in
Applied Earth Observations and Remote Sensing, 8(6):2473–2484, 2015.
D. P. Kingma and M. Welling. Auto-encoding variational bayes. arXiv
preprint arXiv:1312.6114, 2013.
K. Kitada and K. Fukuyama. Land-use and land-cover mapping using a
gradable classification method. Remote Sensing, 4(6):1544–1558, 2012.
V. Kolmogorov. Convergent tree-reweighted message passing for energy
minimization. IEEE transactions on pattern analysis and machine
intelligence, 28(10):1568–1583, 2006.
V. Kolmogorov and R. Zabin. What energy functions can be minimized
via graph cuts? IEEE transactions on pattern analysis and machine
intelligence, 26(2):147–159, 2004.
A. Krizhevsky, I. Sutskever, and G. E. Hinton. Imagenet classification
with deep convolutional neural networks. In Advances in neural
information processing systems, pages 1097–1105, 2012.
J. Li, J. M. Bioucas-Dias, and A.Plaza. Spectral–spatial hyperspectral
image segmentation using subspace multinomial logistic regression and
markov random fields. IEEE Transactions on Geoscience and Remote
Sensing, 50(3):809–823, 2012.
J. Li, J. M. Bioucas-Dias, and A. Plaza. Semisupervised hyperspectral
image segmentation using multinomial logistic regression with active
learning. IEEE Transactions on Geoscience and Remote Sensing,
48(11):4085–4098, 2010.
G. Licciardi, P. R. Marpu, J. Chanussot, and J. A. Benediktsson. Linear
versus nonlinear pca for the classification of hyperspectral data based
on the extended morphological profiles. IEEE Geoscience and Remote
Sensing Letters, 9(3):447–451, 2012.
Z. Lin, Y. Chen, X. Zhao, and G. Wang. Spectral-spatial classification
of hyperspectral image using autoencoders. In IEEE 9th International
Conference on Information, Communications and Signal Processing
(ICICS), 2013, pages 1–5. IEEE, 2013.
L. Ma, M. M. Crawford, and J. Tian. Local manifold learningbased k-nearest-neighbor for hyperspectral image classification. IEEE
Transactions on Geoscience and Remote Sensing, 48(11):4099–4109,
2010.
T. Matsuki, N. Yokoya, and A. Iwasaki. Hyperspectral tree species
classification of japanese complex mixed forest with the aid of lidar
data. IEEE Journal of Selected Topics in Applied Earth Observations
and Remote Sensing, 8(5):2177–2187, 2015.
F. Melgani and L. Bruzzone. Classification of hyperspectral remote
sensing images with support vector machines. IEEE Transactions on
geoscience and remote sensing, 42(8):1778–1790, 2004.
R. Mendez-Rial and J. Martin-Herrero. Efficiency of semi-implicit
schemes for anisotropic diffusion in the hypercube. IEEE Transactions
on Image Processing, 21(5):2389–2398, 2012.
L. Mou, P. Ghamisi, and X. Zhu. Deep recurrent neural networks for
hyperspectral image classification. IEEE Transactions on Geoscience
and Remote Sensing, 2017.
S. Paisitkriangkrai, J. Sherrah, P. Janney, V.-D. Hengel, et al. Effective
semantic pixel labelling with convolutional networks and conditional
random fields. In Proceedings of the IEEE Conference on Computer
Vision and Pattern Recognition Workshops, pages 36–43, 2015.
G. P. Petropoulos, C. Kalaitzidis, and K. P. Vadrevu. Support vector
machines and object-based classification for obtaining land-use/cover
cartography from hyperion hyperspectral imagery. Computers & Geosciences, 41:99–107, 2012.
A. Plaza, J. A. Benediktsson, J. W. Boardman, J. Brazile, L. Bruzzone,
G. Camps-Valls, J. Chanussot, M. Fauvel, P. Gamba, A. Gualtieri,
et al. Recent advances in techniques for hyperspectral image processing.
Remote sensing of environment, 113:S110–S122, 2009.
Y. Qian, M. Ye, and J. Zhou. Hyperspectral image classification based
on structured sparse logistic regression and three-dimensional wavelet
texture features. IEEE Transactions on Geoscience and Remote Sensing,
51(4):2276–2291, 2013.
M. G. Reyes, X. Zhao, D. L. Neuhoff, and T. N. Pappas. Lossy
compression of bilevel images based on markov random fields. In IEEE
International Conference on Image Processing (ICIP), 2007, volume 2,
pages II–373. IEEE, 2007.
D. E. Rumelhart, G. E. Hinton, and R. J. Williams. Learning internal
representations by error propagation. Technical report, DTIC Document,
1985.
H. Z. Shafri, E. Taherzadeh, S. Mansor, and R. Ashurov. Hyperspectral remote sensing of urban areas: an overview of techniques and
applications. Research Journal of Applied Sciences, Engineering and
Technology, 4(11):1557–1565, 2012.
L. Shen and S. Jia. Three-dimensional gabor wavelets for pixel-based
hyperspectral imagery classification. IEEE Transactions on Geoscience
and Remote Sensing, 49(12):5039–5046, 2011.

14

[55] A. Soltani-Farani, H. R. Rabiee, and S. A. Hosseini. Spatial-aware
dictionary learning for hyperspectral image classification. IEEE Transactions on Geoscience and Remote Sensing, 53(1):527–541, 2015.
[56] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov. Dropout: A simple way to prevent neural networks from overfitting.
The Journal of Machine Learning Research, 15(1):1929–1958, 2014.
[57] X. Sun, Q. Qu, N. M. Nasrabadi, and T. D. Tran. Structured priors
for sparse-representation-based hyperspectral image classification. IEEE
Geoscience and Remote Sensing Letters, 11(7):1235–1239, 2014.
[58] C. Szegedy, A. Toshev, and D. Erhan. Deep neural networks for object
detection. In Advances in Neural Information Processing Systems, pages
2553–2561, 2013.
[59] Y. Tarabalka, J. A. Benediktsson, and J. Chanussot. Spectral–spatial
classification of hyperspectral imagery based on partitional clustering
techniques. IEEE Transactions on Geoscience and Remote Sensing,
47(8):2973–2987, 2009.
[60] Y. Tarabalka, M. Fauvel, J. Chanussot, and J. A. Benediktsson. Svm-and
mrf-based method for accurate classification of hyperspectral images.
IEEE Geoscience and Remote Sensing Letters, 7(4):736–740, 2010.
[61] Y. Tarabalka and A. Rana. Graph-cut-based model for spectral-spatial
classification of hyperspectral images. In IEEE International Geoscience
and Remote Sensing Symposium (IGARSS), 2014, pages 3418–3421.
IEEE, 2014.
[62] A. Villa, J. A. Benediktsson, J. Chanussot, and C. Jutten. Hyperspectral
image classification with independent component discriminant analysis.
IEEE transactions on Geoscience and remote sensing, 49(12):4865–
4876, 2011.
[63] J. Xie, L. Xu, and E. Chen. Image denoising and inpainting with deep
neural networks. In Advances in Neural Information Processing Systems,
pages 341–349, 2012.
[64] Y. Xu, Z. Wu, and Z. Wei. Spectral–spatial classification of hyperspectral
image based on low-rank decomposition. IEEE Journal of Selected
Topics in Applied Earth Observations and Remote Sensing, 8(6):2370–
2380, 2015.
[65] J. S. Yedidia, W. T. Freeman, and Y. Weiss. Constructing free-energy
approximations and generalized belief propagation algorithms. IEEE
Transactions on Information Theory, 51(7):2282–2312, 2005.
[66] J. Yue, S. Mao, and M. Li. A deep learning framework for hyperspectral
image classification using spatial pyramid pooling. Remote Sensing
Letters, 7(9):875–884, 2016.
[67] J. Yue, W. Zhao, S. Mao, and H. Liu. Spectral–spatial classification of
hyperspectral images using deep convolutional neural networks. Remote
Sensing Letters, 6(6):468–477, 2015.
[68] H. Zhang, Y. Li, Y. Zhang, and Q. Shen. Spectral-spatial classification
of hyperspectral imagery using a dual-channel convolutional neural
network. Remote Sensing Letters, 8(5):438–447, 2017.
[69] W. Zhu, V. Chayes, A. Tiard, S. Sanchez, D. Dahlberg, A. L. Bertozzi,
S. Osher, D. Zosso, and D. Kuang. Unsupervised classification in
hyperspectral imagery with nonlocal total variation and primal-dual
hybrid gradient algorithm. IEEE Transactions on Geoscience and
Remote Sensing, 55(5):2786–2798, 2017.
[70] D. Zikic, B. Glocker, O. Kutter, M. Groher, N. Komodakis, A. Kamen,
N. Paragios, and N. Navab. Linear intensity-based image registration
by markov random fields and discrete optimization. Medical image
analysis, 14(4):550–562, 2010.

