# Tile-Based Semisupervised Classification of Large-Scale VHR Remote Sensing Images

Authors: Haikel Alhichri; Essam Othman; Mansour Zuair; Nassim Ammour; Yakoub Bazi
Journal: Journal of Sensors
 (2018)

---

## Abstract

This paper deals with the problem of the classification of large-scale very high-resolution (VHR) remote sensing (RS) images in a semisupervised scenario, where we have a limited training set (less than ten training samples per class). Typical pixel-based classification methods are unfeasible for large-scale VHR images. Thus, as a practical and efficient solution, we propose to subdivide the large image into a grid of tiles and then classify the tiles instead of classifying pixels. Our proposed method uses the power of a pretrained convolutional neural network (CNN) to first extract descriptive features from each tile. Next, a neural network classifier (composed of 2 fully connected layers) is trained in a semisupervised fashion and used to classify all remaining tiles in the image. This basically presents a coarse classification of the image, which is sufficient for many RS application. The second contribution deals with the employment of the semisupervised learning to improve the classification accuracy. We present a novel semisupervised approach which exploits both the spectral and spatial relationships embedded in the remaining unlabelled tiles. In particular, we embed a spectral graph Laplacian in the hidden layer of the neural network. In addition, we apply regularization of the output labels using a spatial graph Laplacian and the random Walker algorithm. Experimental results obtained by testing the method on two large-scale images acquired by the IKONOS2 sensor reveal promising capabilities of this method in terms of classification accuracy even with less than ten training samples per class.


---

## Body

## 1. Introduction

The intent of the classification process is to categorize all pixels in a remote sensing (RS) image into one of several land-cover classes. Most classification solutions have used RS images with small sizes (typically less than 1000 × 1000 pixels). However, the latest generation of satellite-based imaging sensors (Pleiades, Sentinel, etc.) acquires big volumes of Earth’s images with high spatial, spectral, and temporal resolution. This leads to images of very large sizes and creates new challenges for land-use classification algorithms. In this case, traditional pixel-based algorithms become unfeasible due to (1) increased computational and memory costs, (2) increased spectral variation within the same land-cover class, and (3) increased variability in the data distributions over this large-scale image, and (4) the need to collect a large set of training samples to properly model the underlying distribution of the data in the image. The latter issue may be caused by the increased resolution which heightens the effect of lighting conditions, the depth dimension, and the sensor’s noise. The same land-covers and even the same objects can be found in images belonging to different classes. Thus, there is a great need to develop efficient solutions to classify large-scale images.One solution to this problem is to use object-based techniques. This can be accomplished first via segmenting the large image into regions and then performing object analysis on the regions. However, segmenting large images is, in itself, still a difficult problem [1]. Another solution is to divide the image into tiles and then perform coarse classification of the image by classifying the tiles instead of pixels as was introduced in [2–4].The problem of tile classification becomes then similar to object classification or recognition in computer vision applications. In those applications, many local descriptors, like bag of words (BOW) [5], local binary patterns (LBP) [6], and histograms of oriented gradients (HOG) [7], with their invariance to geometric and photometric transformations, have been proven effective especially for object recognition. They can be extracted both in sparse and in dense ways. Recently, deep convolutional neural networks (CNN) [8–10] have shown outstanding results in many computer vision applications such as image classification [11], object recognition [12], face recognition [13], medical image analysis [14], speech recognition [15], and traffic flow prediction [16]. Several CNN architectures have been already proposed and tested in computer vision tasks, and most of them have been implemented and made available online [17].CNNs use images directly as input, and they automatically learn a hierarchy of features instead of using handcrafted features. This is accomplished by successively convolving the input image with learned filters to build up a hierarchy of feature maps. The hierarchical approach allows learning more and more complex and highly descriptive set of features. One disadvantages of CNN is the high computational costs during the training phase due to the high number of weights to be estimated (in the range of hundreds of thousands). Another stringent requirement is the need for a huge number of training images [18]. One way to solve this problem is to use CNNs that have been pretrained on huge auxiliary datasets and then transfer the knowledge embedded in them to help extract highly descriptive features in our new domain. These highly descriptive features are then used to train another classifier (such as a smaller neural network) for our particular recognition task.In the remote sensing community, some works have introduced CNN as a solution to classification and object detection problems [19–23]. For the scene classification problem, the reader is referred to this work by Cheng et al. [24] for a recent review of the state-of-the-art review in the area. As a set of sample work, Castelluccio et al. proposed a solution for the classification of RS scenes by starting with a pretrained CNN and then fine-tuning it on the target data in order to improve the accuracy [19]. Zhang et al. proposed a method called gradient boosting random convolutional network (GBRCN), which trains and fuses the result of several CNNs for the purpose of RS scene classification [22]. Chen et al. proposed a method to detect small objects such as vehicles in satellite images using a CNN to extract and learn rich features from the training data [21]. Marmanis et al. extracted feature representations using a pretrained CNN and then fed them into another CNN, which is trained in a supervised manner [20]. Another more recent work proposes a novel new feature for scene classification called bag of convolutional features (BoCF) [23]. This is similar to bag of words concept, except that the codebook is composed of convolutional feature vectors collected from the 5th layer of the VGGNet-16 CNN when it is applied to all training and testing scenes.Another important issue in RS classification is that the collection of a statistically significant and representative amount of training (labelled) samples is a tedious task. Thus, the RS community has introduced semisupervised learning (SSL) methods to tackle this problem [25–29]. SSL methods attempt to exploit the wealth of information provided by unlabelled data, besides the few available labelled data, to improve the performance of the classifier. SSL is based on the assumption that points within the same structure (such as a cluster or a manifold) are likely to have the same label [29]. We can use the large amount of unlabelled data to uncover such structure. For example, in RS images, it is reasonable to assume that samples are likely to have similar labels, if they have close spectral information or if they are neighbours spatially. Many machine-learning researchers have found that unlabelled data, when used in conjunction with a small amount of labelled data, can produce considerable improvement in learning accuracy [30].Semisupervised learning is natural for deep CNN architectures since the latter requires huge training data and in practice, either we never have enough labelled data or it is quite costly to get them. Nonlinear embedding algorithms are popular for use with shallow semisupervised learning techniques such as kernel methods, but they can be applied to deep multilayer architectures as well, as a regularizer either at the output layer or on the hidden layers of the architecture [31]. Graph-based techniques are one way to perform semisupervised learning. In this approach, one tries to preserve the structure present in the unlabelled data by building a graph where each sample data point is a node and the edges encode the similarity between data points. In the literature, graph Laplacian is embedded at the hidden layer and optimized in a joint way within the error function of the network [29]. Another approach is to use the graph Laplacian as a regularizer at the output layer such as in Markov random field (MRF).In this work, we present an effective solution for the classification of large-scale VHR RS images, where the number of labelled samples is limited (no more than ten tiles per class). To deal with the huge size of the VHR images, we subdivide the large image into square tiles and classify them instead of classifying pixels. This is our first contribution in this paper. Our second contribution is presenting a deep architecture for semisupervised classification of large-scale images. This deep architecture is composed of a pretrained CNN, for feature extraction, followed by two fully connected layers for classification. The proposed solution accomplishes semisupervised learning via simultaneous embedding of a graph Laplacian in the hidden layer followed by a spatial graph Laplacian after the output layer for regularization. The experimental results, carried out on two large-scale VHR images, have shown that embedding a graph Laplacian at the feature level and at the output level simultaneously provides significant improvement in the classification result.The rest of the paper is organized as follows: in Section2, we present the details of proposed method. In Section 3, we present the experimental results and discuss them. Finally, we present our concluding remarks and future works in Section 4.


## 2. Materials and Methods

The general overview of our proposed solution is shown in Figure1. The first step in this solution is to divide the large image into a grid of tiles, and the goal is to classify these tiles instead of pixels to produce a coarse tile-based classification of the RS image. Let I be a large VHR RS image that is divided into a grid of Nt tiles Iii=1Nt, where Ii∈Rw×w represents tile i of size w×w. The second step in the algorithm is feature extraction using one of the pretrained CNNs found in the literature. Let xii=1Nt be the set of features extracted from each one of the Nt tiles. Furthermore, let Tr=xi,yii=1n be a training set of labelled sample tiles, where yi∈1,2,…,C are the corresponding class labels.Figure 1
Overview of the proposed deep CNN solution for the semisupervised classification of large RS images.The classification step is carried out using a fully connected neural network (NN), as shown in Figure1, with one hidden layer and an output (softmax) layer. As mentioned previously, using this proposed architecture, we can implement semisupervised learning in different ways. We can embed the graph Laplacian at the hidden layer of the fully connected neural network (i.e., at the feature level). We can also embed a graph Laplacian at the output level, which implements what is called output regularization. Finally, we can implement a suitable combination thereof. Our final proposed solution shown in Figure 1 embeds a graph Laplacian (spectral or spatial) in the hidden layer of the NN and also adds an additional layer that performs regularization of the output labels. In the remainder of this section, we explain how to use a pretrained CNN for feature extraction. Then the next subsection shows how to perform classification using a fully connected neural network. Finally, we explain how to achieve semisupervised learning via embedding a graph Laplacian at the hidden layer of the neural network and also embedding another graph Laplacian to perform output regularization spatially.


### 2.1. Feature Extraction from Tiles

The second step in the proposed solution is to extract feature descriptors from the image tiles. Several types of handcrafted feature descriptors are presented in the literature for image representation, such as bag of words (BOW) [5], local binary patterns (LBP) [6], or histograms of oriented gradients (HOG) [7]. However, more recently, the new state-of-the-art approach is to use feature descriptors that can be learned instead of handcrafted features. CNNs play a major role in this regard as they have been proven effective in feature learning from images. However, as CNN needs a huge amount of training data, one can resort to pretrained CNN as a compromise solution. In this case, a pretrained CNN becomes like a feature extraction tool instead of a classifier as it is originally designed. This approach exploits the huge information embedded in pretrained CNNs to provide highly descriptive features from the tiles.A CNN is composed of several layers that do different functions such as convolution, pooling, activation function, and classification as shown in Figure1. The last layer is always a softmax layer which is used for classification. Given this deep nature, one can extract features at different layers of the CNN. In this work, we take the output of the hidden layer before the softmax layer as feature vector representation. Thus, given a CNN with L layers, we feed each tile image Ii, i=1,…,Nt as input to the CNN and generate a feature representation xi∈RD at layer k, of the CNN, with k=L−1, that is,
(1)xi=fkCNN…f2CNNf1CNNIi,i=1,…,Nt,where D is the dimensionality of the feature vector.


### 2.2. Classification via a Fully Connected Layer

Given a training setTr=xi,yii=1n composed of feature vectors and their corresponding labels, the hidden layer takes the input xi and maps it to another representation hi1∈Rd1 of dimension d1 through the nonlinear activation function f as follows:
(2)hi1=fW1xi,where W1∈Rd1×d is the mapping weight matrix. A typical choice of the activation function is the sigmoid function, that is, fv=1/1+exp−v. For simplicity, we omit the bias vector in the expression as it can be incorporated as an additional column vector in the mapping matrix; then in that case, the feature vector should be appended by the value 1.Thesoftmax layer performs multiclass classification and takes as input the resulting hidden representation hi1. It produces, as a result, an estimate of the posterior probability for each class label c=1,2,…,C as follows:
(3)py^i=c∣xi=expwc2Thi1∑j=1Cexpwj2Thi1,where W2=w12w22…wC2∈Rd×C are the weights of the softmax layer and the superscript ∙T refers to the transpose operation. To learn the weights θ=W1,W2 representing the complete network structure, we minimize the training error on the labelled training data. The cost function is then formulated as follows:
(4)Jθ,D1=argminEnetD1,where Dl is the set of labelled data and EnetDl is the cross-entropy loss, which measures the error between the actual network outputs and the desired outputs of the labelled source data. As the outputs of the network are probabilistic, we propose to maximize the log-posterior probability to learn the network weights, which is equivalent to minimizing the so-called cross-entropy error:
(5)EnetDl=−1n∑i=1ns∑k=1K1yi=klnexpwk2Thk1∑j=1Kexpwj2Thl1,where 1∙ is an indicator function that takes 1 if the statement is true otherwise it takes 0 and the superscript T refers to matrix transpose.


### 2.3. Graph-Based Semisupervised Learning

The idea of graph-based semisupervised learning is to build a graph that connects similar sample data points to each other, and then use this graph to propagate estimated labels among similar data points. This follows directly from the assumption that similar data points should have similar labels, which the graph tries to encode.Figure2 illustrates how a Laplacian matrix is built using a simple example. Let G=V,E be a graph with vertices v∈V and edges e∈E. Let eij denote an edge spanning two vertices vi and vj and ωij is the associated weight. The degree of a vertex vi is di=∑ωij for all edges eij incident on vi. Then the tile-based Laplacian matrix indexed by vertices vi and vj is given by:
(6)Lij=dij,ifi=j,−ωij,ifiandjareconnected,0,otherwise.Figure 2
Illustration of graph Laplacian building. (a) Sample graph with 6 vertices and 7 edges. (b) Laplacian matrix. (c) Spectral graph Laplacian. Each tile (a node in the graph) is connected to thek-most similar tiles. (d) Spatial graph Laplacian, each tile (a node in the graph) is connected to its 8 neighbours.In our case, each tile image is associated with a vertex in the graph, while edges are defined in several ways. Thek-nearest graph, for example, connects each vertex to the k-most similar vertices. The number k is usually left as a parameter to be determined. Another approach is to connect each vertex to all vertices whose similarity is below a certain threshold value. This is known as a spectral graph and is illustrated in Figure 2(d).However, here, the tiles are also spatially dependent. Thus, we can also consider adjacency as a similarity, that is, we can make the assumption that neighbouring tiles are similar and should have similar labels. In this case, the number of connections per node is usually limited to two choices: 8 or 4 depending on whether we want to consider all 8 neighbouring tiles or only 4. This is known as a spatial graph and is illustrated in Figure2(c).Finally, a weightwij is associated with each edge eij in the graph. This weight is defined as the similarity between image tiles corresponding to vertices vi and vj. It can be calculated using many types of similarity measures or distance measures. A common choice for obtaining these weights is the Gaussian weighting function, that is, ωjk=exp−βxj−xk2 and β is a free parameter (usually set to 1).


#### 2.3.1. Embedding Graph Laplacian in the Hidden Layer

The first approach to implement semisupervised learning is by embedding a spectral or spatial graph Laplacian in the hidden layer. The fully connected neural network is supplied with labelled and unlabelled data. Then to learn the weightsθ=W1,W2 representing the complete network structure, we propose to simultaneously minimize: (i) the training error on the labelled data and (ii) the energy of the graph Laplacian built on all unlabelled tiles. The proposed cost function is then formulated as follows:
(7)Jθ,D1,Du=arg minEnetDl+λ1LapDu,where λ1 is a regularization parameter. The first term is the same as (4), while the second term is added for graph regularization. This basically minimizes the energy of the graph and can be written as follows:
(8)LapDu=12∑ejk∈EωjkW1xj−W1xk22.It can be shown that (8) can be written in the following matrix form:
(9)LapDu=12TraceW1XuLW1XuT,where Xu=x1⋯xnu∈Rnu×d1 represents the unlabelled samples.By substituting the terms in (5) and (9), the total cost function Jθ is then given by
(10)Jθ,D1,Du=arg min−1n∑i=1nl∑k=1K1yi=klnexpwk2Thk1∑l=1Kexpwl2Thl1+λ12TraceW1XuLW1XuT.


#### 2.3.2. Spatial Regularization of the Output Layer

Another approach to exploit the graph Laplacian is an output regularizer based on the random Walker (RW) algorithm. This technique is similar to a Markov random field (MRF) approach; however, unlike the MRF, it yields a closed-form solution as opposed to an iterative one [32]. This regularization step needs a spatial graph Laplacian and posterior probability estimates as initial inputs. Then the RW algorithm can improve the classification result by minimizing an energy function of the graph Laplacian, defined as follows:
(11)minf∑eij∈Eωijfi−fj2+λ2∑vi∈Vμifi−fi∗2,where fi and fj are now probability vectors (of dimension C) associated with vertices (tiles) vi and vj, respectively. The first term is related to image smoothness since it computes the norm between pixels. The second term is related to data fidelity since it computes a norm between the initial estimation fi∗ (provided by the fully connected NN) and the new label estimatesfi. Finally, μi and λ2 are local and global weights, respectively, enforcing that fidelity. Moreover, the local weights μi could be set equal to di, but for simplicity, we suppose μi=1 for all tiles in the image.The set of vertices of the graph can be written asV=VT∪VU, where VT and VU are the set of labelled and unlabelled vertices, respectively. In this case, the Laplacian matrix L can be written as L=LTBBLU. After differentiating, setting the result to 0 and solving the matrix equation, we found that the probabilities associated with the unlabelled tiles are obtained by the following closed-form solution:
(12)FU=LU+λ2I−1−BTFT+λFU∗,where
(13)F=FTFUand FU∗ are the initial estimated probabilities (posterior probabilities of the fully connected NN classifier) associated with the unlabelled tiles. Note that LU+λ2I is a positive semidefinite matrix, and L is a sparse (Laplacian) matrix. Thus, the linear system of equations in (12) can be solved in linear time.


## 2.1. Feature Extraction from Tiles

The second step in the proposed solution is to extract feature descriptors from the image tiles. Several types of handcrafted feature descriptors are presented in the literature for image representation, such as bag of words (BOW) [5], local binary patterns (LBP) [6], or histograms of oriented gradients (HOG) [7]. However, more recently, the new state-of-the-art approach is to use feature descriptors that can be learned instead of handcrafted features. CNNs play a major role in this regard as they have been proven effective in feature learning from images. However, as CNN needs a huge amount of training data, one can resort to pretrained CNN as a compromise solution. In this case, a pretrained CNN becomes like a feature extraction tool instead of a classifier as it is originally designed. This approach exploits the huge information embedded in pretrained CNNs to provide highly descriptive features from the tiles.A CNN is composed of several layers that do different functions such as convolution, pooling, activation function, and classification as shown in Figure1. The last layer is always a softmax layer which is used for classification. Given this deep nature, one can extract features at different layers of the CNN. In this work, we take the output of the hidden layer before the softmax layer as feature vector representation. Thus, given a CNN with L layers, we feed each tile image Ii, i=1,…,Nt as input to the CNN and generate a feature representation xi∈RD at layer k, of the CNN, with k=L−1, that is,
(1)xi=fkCNN…f2CNNf1CNNIi,i=1,…,Nt,where D is the dimensionality of the feature vector.


## 2.2. Classification via a Fully Connected Layer

Given a training setTr=xi,yii=1n composed of feature vectors and their corresponding labels, the hidden layer takes the input xi and maps it to another representation hi1∈Rd1 of dimension d1 through the nonlinear activation function f as follows:
(2)hi1=fW1xi,where W1∈Rd1×d is the mapping weight matrix. A typical choice of the activation function is the sigmoid function, that is, fv=1/1+exp−v. For simplicity, we omit the bias vector in the expression as it can be incorporated as an additional column vector in the mapping matrix; then in that case, the feature vector should be appended by the value 1.Thesoftmax layer performs multiclass classification and takes as input the resulting hidden representation hi1. It produces, as a result, an estimate of the posterior probability for each class label c=1,2,…,C as follows:
(3)py^i=c∣xi=expwc2Thi1∑j=1Cexpwj2Thi1,where W2=w12w22…wC2∈Rd×C are the weights of the softmax layer and the superscript ∙T refers to the transpose operation. To learn the weights θ=W1,W2 representing the complete network structure, we minimize the training error on the labelled training data. The cost function is then formulated as follows:
(4)Jθ,D1=argminEnetD1,where Dl is the set of labelled data and EnetDl is the cross-entropy loss, which measures the error between the actual network outputs and the desired outputs of the labelled source data. As the outputs of the network are probabilistic, we propose to maximize the log-posterior probability to learn the network weights, which is equivalent to minimizing the so-called cross-entropy error:
(5)EnetDl=−1n∑i=1ns∑k=1K1yi=klnexpwk2Thk1∑j=1Kexpwj2Thl1,where 1∙ is an indicator function that takes 1 if the statement is true otherwise it takes 0 and the superscript T refers to matrix transpose.


## 2.3. Graph-Based Semisupervised Learning

The idea of graph-based semisupervised learning is to build a graph that connects similar sample data points to each other, and then use this graph to propagate estimated labels among similar data points. This follows directly from the assumption that similar data points should have similar labels, which the graph tries to encode.Figure2 illustrates how a Laplacian matrix is built using a simple example. Let G=V,E be a graph with vertices v∈V and edges e∈E. Let eij denote an edge spanning two vertices vi and vj and ωij is the associated weight. The degree of a vertex vi is di=∑ωij for all edges eij incident on vi. Then the tile-based Laplacian matrix indexed by vertices vi and vj is given by:
(6)Lij=dij,ifi=j,−ωij,ifiandjareconnected,0,otherwise.Figure 2
Illustration of graph Laplacian building. (a) Sample graph with 6 vertices and 7 edges. (b) Laplacian matrix. (c) Spectral graph Laplacian. Each tile (a node in the graph) is connected to thek-most similar tiles. (d) Spatial graph Laplacian, each tile (a node in the graph) is connected to its 8 neighbours.In our case, each tile image is associated with a vertex in the graph, while edges are defined in several ways. Thek-nearest graph, for example, connects each vertex to the k-most similar vertices. The number k is usually left as a parameter to be determined. Another approach is to connect each vertex to all vertices whose similarity is below a certain threshold value. This is known as a spectral graph and is illustrated in Figure 2(d).However, here, the tiles are also spatially dependent. Thus, we can also consider adjacency as a similarity, that is, we can make the assumption that neighbouring tiles are similar and should have similar labels. In this case, the number of connections per node is usually limited to two choices: 8 or 4 depending on whether we want to consider all 8 neighbouring tiles or only 4. This is known as a spatial graph and is illustrated in Figure2(c).Finally, a weightwij is associated with each edge eij in the graph. This weight is defined as the similarity between image tiles corresponding to vertices vi and vj. It can be calculated using many types of similarity measures or distance measures. A common choice for obtaining these weights is the Gaussian weighting function, that is, ωjk=exp−βxj−xk2 and β is a free parameter (usually set to 1).


### 2.3.1. Embedding Graph Laplacian in the Hidden Layer

The first approach to implement semisupervised learning is by embedding a spectral or spatial graph Laplacian in the hidden layer. The fully connected neural network is supplied with labelled and unlabelled data. Then to learn the weightsθ=W1,W2 representing the complete network structure, we propose to simultaneously minimize: (i) the training error on the labelled data and (ii) the energy of the graph Laplacian built on all unlabelled tiles. The proposed cost function is then formulated as follows:
(7)Jθ,D1,Du=arg minEnetDl+λ1LapDu,where λ1 is a regularization parameter. The first term is the same as (4), while the second term is added for graph regularization. This basically minimizes the energy of the graph and can be written as follows:
(8)LapDu=12∑ejk∈EωjkW1xj−W1xk22.It can be shown that (8) can be written in the following matrix form:
(9)LapDu=12TraceW1XuLW1XuT,where Xu=x1⋯xnu∈Rnu×d1 represents the unlabelled samples.By substituting the terms in (5) and (9), the total cost function Jθ is then given by
(10)Jθ,D1,Du=arg min−1n∑i=1nl∑k=1K1yi=klnexpwk2Thk1∑l=1Kexpwl2Thl1+λ12TraceW1XuLW1XuT.


### 2.3.2. Spatial Regularization of the Output Layer

Another approach to exploit the graph Laplacian is an output regularizer based on the random Walker (RW) algorithm. This technique is similar to a Markov random field (MRF) approach; however, unlike the MRF, it yields a closed-form solution as opposed to an iterative one [32]. This regularization step needs a spatial graph Laplacian and posterior probability estimates as initial inputs. Then the RW algorithm can improve the classification result by minimizing an energy function of the graph Laplacian, defined as follows:
(11)minf∑eij∈Eωijfi−fj2+λ2∑vi∈Vμifi−fi∗2,where fi and fj are now probability vectors (of dimension C) associated with vertices (tiles) vi and vj, respectively. The first term is related to image smoothness since it computes the norm between pixels. The second term is related to data fidelity since it computes a norm between the initial estimation fi∗ (provided by the fully connected NN) and the new label estimatesfi. Finally, μi and λ2 are local and global weights, respectively, enforcing that fidelity. Moreover, the local weights μi could be set equal to di, but for simplicity, we suppose μi=1 for all tiles in the image.The set of vertices of the graph can be written asV=VT∪VU, where VT and VU are the set of labelled and unlabelled vertices, respectively. In this case, the Laplacian matrix L can be written as L=LTBBLU. After differentiating, setting the result to 0 and solving the matrix equation, we found that the probabilities associated with the unlabelled tiles are obtained by the following closed-form solution:
(12)FU=LU+λ2I−1−BTFT+λFU∗,where
(13)F=FTFUand FU∗ are the initial estimated probabilities (posterior probabilities of the fully connected NN classifier) associated with the unlabelled tiles. Note that LU+λ2I is a positive semidefinite matrix, and L is a sparse (Laplacian) matrix. Thus, the linear system of equations in (12) can be solved in linear time.


## 2.3.1. Embedding Graph Laplacian in the Hidden Layer

The first approach to implement semisupervised learning is by embedding a spectral or spatial graph Laplacian in the hidden layer. The fully connected neural network is supplied with labelled and unlabelled data. Then to learn the weightsθ=W1,W2 representing the complete network structure, we propose to simultaneously minimize: (i) the training error on the labelled data and (ii) the energy of the graph Laplacian built on all unlabelled tiles. The proposed cost function is then formulated as follows:
(7)Jθ,D1,Du=arg minEnetDl+λ1LapDu,where λ1 is a regularization parameter. The first term is the same as (4), while the second term is added for graph regularization. This basically minimizes the energy of the graph and can be written as follows:
(8)LapDu=12∑ejk∈EωjkW1xj−W1xk22.It can be shown that (8) can be written in the following matrix form:
(9)LapDu=12TraceW1XuLW1XuT,where Xu=x1⋯xnu∈Rnu×d1 represents the unlabelled samples.By substituting the terms in (5) and (9), the total cost function Jθ is then given by
(10)Jθ,D1,Du=arg min−1n∑i=1nl∑k=1K1yi=klnexpwk2Thk1∑l=1Kexpwl2Thl1+λ12TraceW1XuLW1XuT.


## 2.3.2. Spatial Regularization of the Output Layer

Another approach to exploit the graph Laplacian is an output regularizer based on the random Walker (RW) algorithm. This technique is similar to a Markov random field (MRF) approach; however, unlike the MRF, it yields a closed-form solution as opposed to an iterative one [32]. This regularization step needs a spatial graph Laplacian and posterior probability estimates as initial inputs. Then the RW algorithm can improve the classification result by minimizing an energy function of the graph Laplacian, defined as follows:
(11)minf∑eij∈Eωijfi−fj2+λ2∑vi∈Vμifi−fi∗2,where fi and fj are now probability vectors (of dimension C) associated with vertices (tiles) vi and vj, respectively. The first term is related to image smoothness since it computes the norm between pixels. The second term is related to data fidelity since it computes a norm between the initial estimation fi∗ (provided by the fully connected NN) and the new label estimatesfi. Finally, μi and λ2 are local and global weights, respectively, enforcing that fidelity. Moreover, the local weights μi could be set equal to di, but for simplicity, we suppose μi=1 for all tiles in the image.The set of vertices of the graph can be written asV=VT∪VU, where VT and VU are the set of labelled and unlabelled vertices, respectively. In this case, the Laplacian matrix L can be written as L=LTBBLU. After differentiating, setting the result to 0 and solving the matrix equation, we found that the probabilities associated with the unlabelled tiles are obtained by the following closed-form solution:
(12)FU=LU+λ2I−1−BTFT+λFU∗,where
(13)F=FTFUand FU∗ are the initial estimated probabilities (posterior probabilities of the fully connected NN classifier) associated with the unlabelled tiles. Note that LU+λ2I is a positive semidefinite matrix, and L is a sparse (Laplacian) matrix. Thus, the linear system of equations in (12) can be solved in linear time.


## 3. Results

### 3.1. Dataset Description

We tested the proposed solution on two large-scale VHR RS images. The first image is of size 6500 × 10400 taken over the city of Jeddah in Saudi Arabia by the IKONOS-2 sensor in July 2004 (see Figure3(a)). The image has three spectral bands with a spatial resolution of 1 m. It has been divided into 3977 tiles of size 128 × 128, of which 1949 tiles are unlabelled, while 2028 are labelled into nine land-cover types including asphalt (73 tiles), bare soil (61 tiles), beach (60 tiles), cargo (174 tiles), dense residential (203 tiles), medium residential (193 tiles), sparse residential (47 tiles), trees (74 tiles), and water (1216 tiles).Figure 3
Large VHR dataset used and its qualitative results. (a) Original image of size 6500 × 10400 with example training tiles (10 per class), (b) testing map (ground truth), (c) legend, (d) NN classification result (OA = 89%), (e) final map (OA = 92.2%), and (f) tiles corrected by semisupervised learning.
(a)
(b)
(c)
(d)
(e)
(f)The second image is of size 14000 × 14000 and is taken from the city of Riyadh by GeoEye sensor in June 2010 (see Figure4(a)). The image has three spectral bands with a spatial resolution of 0.5 m. It has been divided into 4900 tiles of which 3860 are unlabelled, while 1040 are labelled into five land-cover types including bare soil (163 tiles), freeway (185 tiles), parking (109 tiles), residential (486 tiles), and sparse residential (97 tiles).Figure 4
Large VHR dataset used and its qualitative results. (a) Original image of size 13440 × 13440 with example training tiles (10 per class), (b) testing map (ground truth), (c) legend, (d) NN classification result (OA = 91%), (e) final map (OA = 96.5%), and (f) tiles corrected by semisupervised learning.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.2. Experimental Setup

Both large images have been subdivided into a grid of equal size tiles. The Jeddah image is divided in tiles of size 128 × 128 pixels, while the Riyadh image tiles are of size 196 × 196 pixels. We start the division from the top left corner, and any leftover pixels (in case the size of the large image is not a multiple of 128) appear at the right and bottom sides of the image. The ground truth classification maps are created manually by inspecting a random set of tiles (the test set) one by one and assigning one of the 9 labels to all tiles in the test set. Testing maps for both the Jeddah and Riyadh images are shown in Figures5(b) and 6(b), respectively. For both images, we extract feature representations from all tiles including HOG, LBP, BOW, and several pretrained CNN architectures. The set of pretrained CNNs used includes the imagenet-vgg-very-deep-16 and imagenet-vgg-m models by the VGG group at the University of Oxford [33, 34], the GoogleNet model developed by a team from Google Inc. [35], and Microsoft’s ResNet-152 model [36]. These CNNs are selected because they have won the first and second places in the ImageNet Large-Scale Visual Recognition Challenge- (ILSVRC-) 2014 and 2015 challenges. Table 1 shows the sizes of each feature representation.Figure 5
Results of output regularization via spatial graph Laplacian with sensitivity analysis of the parameterλ2 in (11), (a) for Riyadh dataset and (b) for the Jeddah dataset.
(a)
(b)Figure 6
Results of embedding a graph Laplacian in the hidden layer with sensitivity analysis of the parameterλ1 in (10). (a) Riyadh dataset with spectral graph. (b) Jeddah dataset with spectral graph. (c) Riyadh dataset with spatial graph. (d) Jeddah dataset with spatial graph.
(a)
(b)
(c)
(d)Table 1
Parameters used by various feature extraction techniques.
Feature type
Feature vector length
Note
HOG
4365
LBP
256
8 neighbourhood, uniform pattern
BOW
4000
Vocabulary size = 4000
CNN—vgg-very-deep-16
4096
CNN—GoogLeNet
1024
CNN—vgg-m
1024
CNN—ResNet-152
2048As for the assessment method, we use the Overall Accuracy (OA) and Average Accuracy (AA) measures. To assess the statistical significance of the improvements achieved by applying the semisupervised learning, we computed McNemar’s test, which is based on the standardized normal test statistic:(14)Zij=fij−fjifij−fji,where Zij measures the pairwise statistical significance of the difference between the accuracies of the ith and jth classifiers. The variable fij stands for the number of samples classified correctly and wrongly by the ith and jth classifier, respectively. Thus, fij and fji are the counts of the classified samples on which the two classifiers disagree. At the commonly used 5% level of significance, the improvement of classifier j over i is said to be significant if Zij<−1.96. Finally, all experiments are repeated 10 times with different random training tiles and then the average accuracy values are presented.


### 3.3. Comparison of Various Pretrained CNN and Finding the Best Hidden Layer Size

In this first experiment, we compare the performance of feature representations obtained by the pretrained CNNs and the ones obtained by BOW, LBP, and HOG. We also investigate in this experiment the sensitivity of the results with respect to the hidden layer size for all the types of feature descriptors considered in this study. We train a fully connected NN with the following hidden layer sizes 32, 64, 128, 256, 512, and 1024. We perform the experiment 10 times; each time, we select a different set of training sample tiles per class randomly and average the accuracy results.Figure7 shows the overall accuracies (OA) achieved for each type of feature descriptor and with different hidden layer sizes. It is clear from this experiment that CNN features are quite superior to all other types of handcrafted features regardless of the hidden layer size. This confirms other results about CNN that is reported in the object recognition literature.Figure 7
Comparison between various feature descriptors in terms of overall accuracy (OA). (a) OA for the Jeddah-IKONOS2 dataset and (b) OA for the Riyadh-GEOEYE dataset.
(a)
(b)Furthermore, if we consider CNN type features only, VGG very-deep-16 pretrained CNN, in particular, provided the best OA. We can also conclude that in general the range [64, 256] constitutes the best values for the hidden layer size. For the remainder of this paper, the hidden layer size of 256 and the VGG very-deep-16 pretrained CNN will be used in all experiments.


### 3.4. Semisupervised Learning Results

In this set of experiments, we present the results of applying semisupervised learning via graph Laplacian. We first implement embedding a spatial graph Laplacian as an output regularizer based on the RW algorithm without embedding any graph in the hidden layer. This step has a regularization parameterλ2 as shown in (11). In Figure 5, we show the Z score achieved for different values of the regularization parameter in the range of 1…20. From these results, we learn that output regularization does not provide significant improvement for the Riyadh dataset as the Z score is greater than −1.96.However, the improvement for the Jeddah dataset is significant (<−1.96). We also learn that the best value for the parameterλ2 is around 4 to 6; thus, we decided to fix this parameter to 5 for optimal performance of the method.In the second experiment, we present the results of embedding a graph Laplacian in the hidden layer. We investigated embedding both a spectral graph Laplacian and spatial graph Laplacian. The cost function to be optimized, shown in (10), has a regularization parameterλ1. Figure 6 shows the Z score achieved for different values of the regularization parameter in the range of 0.04…0.15. From these results, we learn that embedding a spatial graph Laplacian in the hidden layer degrades the classification accuracy significantly as the Z scores in Figures 6(c) and 6(d) are all above −1.96. On the other hand, embedding a spectral graph Laplacian in the hidden layer does provide significant improvement for both datasets as the Z score is less than −1.96. As for the best value for the parameterλ1, it is clear, from Figures 6(a) and 6(b), that it is around the values 0.08 to 0.1; thus, we decided to fix this parameter to 0.1 for optimal performance of the method.Next, we present the results of the complete system shown in Figure1. The results for the Jeddah and Riyadh image are shown in Tables 2 and 3, respectively. These tables show the OA, AA, and per class accuracy for 5 different scenarios: (1) using the fully connected NN only, (2) using the NN followed by spatial regularization using RW algorithm, (3) using the NN with spatial graph Laplacian in the hidden layer, (4) using the NN with spectral graph Laplacian in the hidden layer, and finally, (5) using the NN with both spectral graph Laplacian in the hidden layer followed by spatial regularization using RW algorithm, which is the complete solution proposed in this paper. These tables also show execution times in seconds, which include both training and testing times.Table 2
Classification results for the Jeddah dataset using one hidden layer of size 256.
OA
AA
Time (s)
Z score
Class 1
Class 2
Class 3
Class 4
Class 5
Class 6
Class 7
Class 8
Class 9
VGG-m
NN
86.05
72.84
0.23
0.00
60.32
82.35
54.40
82.38
56.58
73.77
77.03
71.72
96.97
NN + RW
86.42
73.24
0.04
−1.23
60.32
81.76
53.60
84.45
53.68
76.56
77.57
73.91
97.28
Spatial graph in the hidden layer
58.24
36.42
19.86
20.96
17.30
61.57
18.40
8.36
21.86
52.29
36.49
34.22
77.33
Spectral graph in the hidden layer
87.75
77.41
25.29
−2.71
80.00
84.90
57.40
81.77
80.57
59.34
86.22
69.69
96.81
NN + spectral graph + RW
88.33
78.08
0.04
−3.59
80.64
84.71
56.80
83.05
81.25
61.97
85.95
71.41
97.01
GoogleNet
NN
84.48
73.23
0.23
0.00
71.75
87.45
55.20
71.89
75.86
44.81
83.78
72.35
96.00
NN + RW
85.21
74.39
0.04
−2.44
72.86
87.45
55.20
75.80
76.94
46.12
84.59
74.38
96.14
Spatial graph in the hidden layer
77.37
58.45
20.05
6.59
47.30
87.65
11.60
78.72
70.21
35.57
58.92
45.31
90.81
Spectral graph in the hidden layer
85.97
76.82
26.38
−2.57
73.65
91.57
60.60
81.10
87.82
39.07
91.08
70.94
95.55
NN + spectral graph + RW
86.46
77.76
0.04
−3.32
74.92
91.77
60.80
84.27
89.33
37.60
91.89
73.59
95.67
ResNet
NN
87.68
79.33
0.37
0.00
81.59
83.14
63.80
91.04
75.13
70.00
81.35
73.59
94.36
NN + RW
88.33
79.69
0.05
−2.59
81.27
84.51
62.60
91.28
75.49
70.22
82.70
73.91
95.26
Spatial graph in the hidden layer
78.82
62.59
33.94
10.27
46.03
93.34
33.80
82.38
61.40
64.86
39.73
52.97
88.78
Spectral graph in the hidden layer
88.65
82.96
38.38
−2.00
84.76
91.96
78.00
91.95
79.22
69.73
86.76
70.16
94.12
NN + spectral graph + RW
89.01
82.92
0.06
−2.63
84.29
92.16
76.60
92.44
79.64
69.67
87.30
69.53
94.70
Very-deep-16
NN
88.41
78.52
0.62
0.00
74.76
88.63
67.60
96.34
57.67
76.45
82.16
65.63
97.45
NN + RW
89.37
78.98
0.09
−3.77
80.95
88.24
84.00
93.90
50.26
86.89
81.08
51.56
96.85
Spatial graph in the hidden layer
68.79
56.74
60.98
16.15
44.60
67.06
29.60
54.15
67.41
48.64
72.16
49.06
77.97
Spectral graph in the hidden layer
90.28
83.45
96.01
−3.00
82.86
90.59
74.80
93.96
74.77
73.28
95.67
68.28
96.86
NN + spectral graph + RW
91.23
83.68
0.09
−4.67
84.13
90.59
72.40
94.39
74.51
74.48
95.13
69.22
98.25Table 3
Classification results for the Riyadh dataset using one hidden layer of size 256.
OA
AA
Time (s)
Z score
Class 1
Class 2
Class 3
Class 4
Class 5
VGG-m
NN
87.96
81.98
0.22
—
99.09
71.49
86.26
96.53
56.55
NN + RW
88.07
82.21
0.05
−0.49
99.15
71.37
87.27
96.56
56.67
Spatial graph in the hidden layer
81.52
71.48
24.00
6.25
81.90
74.57
78.79
94.33
27.82
Spectral graph in the hidden layer
93.12
88.98
30.36
−5.30
99.22
85.03
91.01
98.61
71.04
NN + spectral graph + RW
93.18
89.07
0.05
−5.38
99.22
85.09
91.52
98.64
70.92
GoogleNet
NN
89.03
88.63
0.37
—
94.51
85.89
82.22
89.48
91.03
NN + RW
89.43
89.05
0.08
−1.06
94.51
85.43
83.64
90.09
91.61
Spatial graph in the hidden layer
80.97
75.54
23.97
5.76
96.21
77.77
70.61
85.50
47.59
Spectral graph in the hidden layer
93.45
91.03
55.17
−4.80
97.91
90.34
87.98
96.26
82.64
NN + spectral graph + RW
93.84
91.40
0.09
−5.20
97.91
90.40
88.18
96.83
83.68
ResNet
NN
89.58
86.33
0.35
—
88.30
81.09
87.78
95.40
79.08
NN + RW
89.72
86.44
0.07
−0.50
89.22
81.26
86.97
95.46
79.31
Spatial graph in the hidden layer
73.64
65.76
39.26
8.61
99.68
30.91
90.20
87.19
20.81
Spectral graph in the hidden layer
95.13
93.56
45.94
−5.86
98.76
89.54
96.36
97.50
85.63
NN + spectral graph + RW
95.29
93.80
0.07
−6.00
99.09
89.43
96.47
97.61
86.44
Very-deep-16
NN
90.97
88.09
0.63
—
95.36
75.60
94.85
96.81
77.82
NN + RW
93.13
91.84
0.11
−1.71
97.39
82.29
95.96
96.22
87.36
Spatial graph in the hidden layer
79.73
67.40
69.42
6.74
80.65
62.29
83.64
97.19
13.22
Spectral graph in the hidden layer
94.72
92.21
91.27
−4.72
98.00
88.85
95.12
98.29
80.79
NN + spectral graph + RW
95.26
92.91
0.11
−5.40
98.50
90.06
94.45
98.55
82.99The results show clearly that semisupervised learning has improved the accuracy significantly for all types of CNN. However, in agreement with an earlier experiment, the best results achieved are with thevgg-very-deep-16 pretrained CNN, with an OA of 91.23% and an AA of 83.68%. Also, the semisupervised learning has significantly improved the accuracy of the NN as the Z score reached −4.67. Similarly for the Riyadh image, the results shown in Table 3 indicate an even better improvement induced by semisupervised learning with a Z score below −5 for all types of CNN. Again, the best results are achieved using the vgg-very-deep-16 pretrained CNN, with an OA of 95.26%, an AA of 92.91%, and a Z score of −5.40.Finally, qualitative results are also shown in Figures3 and 4 for the Jeddah and Riyadh datasets, respectively. Here, we show the classification maps that have achieved the highest OA after semisupervised learning among the 10 successive runs of the algorithm. In Figures 3(e) and 4(e), we show the tiles whose labels are corrected by semisupervised learning. Comparing Figure 3(d) to Figure 3(c), one can see the type of improvement produced by semisupervised learning, namely, the reduced noise particularly in the water and in residential classes. As for the Riyadh classification maps shown in Figures 4(c) and 4(d), the improvements are clearly visible in the freeway class, where it is more recognizable as roads.In this final experiment, we compare the proposed method to a state-of-the-art method in the literature called the Laplacian extreme learning machine (LapELM) method [37]. To make the comparison fairer, we have changed the features used in the original method [37] to the same CNN features used here. Thus, the two methods now use the same CNN features extracted with the help of the vgg-very-deep-16 pretrained CNN. As for the other parameters of the LapELM method, we have set the ranges for the regularization parameter C to [0.001, 1000] and the RBF kernel parameter γ to [0.01, 5]. To estimate the best parameter values in these ranges, we use the deferential evolution technique, which outperforms the grid search technique as shown in [37]. As for the proposed semisupervised method, we note that the batch size parameter is usually set to 100 in the literature. However, here, we should take care not to set the batch size to be larger than the total number of training samples (number of samples per class × the number of classes). For this experiment, we set it equal to the number of samples per class times the number of classes.We also study the sensitivity of the method to reduced training samples (less than 10). Figure8 shows the OA for sample sizes from 2 to 10 for both methods being compared. Figure 8(a) shows the results for the Jeddah dataset, while Figure 8(b) shows the result for the Riyadh dataset. The first observation is that the classification accuracy for both methods degrades gracefully when the number of training samples per class decreases. However, the results confirm again the superiority of the proposed method compared to Laplacian ELM.Figure 8
Comparison between the proposed semisupervised NN and the Laplacian ELM method. (a) Results for Jeddah dataset. (b) Results for Riyadh dataset.
(a)
(b)


## 3.1. Dataset Description

We tested the proposed solution on two large-scale VHR RS images. The first image is of size 6500 × 10400 taken over the city of Jeddah in Saudi Arabia by the IKONOS-2 sensor in July 2004 (see Figure3(a)). The image has three spectral bands with a spatial resolution of 1 m. It has been divided into 3977 tiles of size 128 × 128, of which 1949 tiles are unlabelled, while 2028 are labelled into nine land-cover types including asphalt (73 tiles), bare soil (61 tiles), beach (60 tiles), cargo (174 tiles), dense residential (203 tiles), medium residential (193 tiles), sparse residential (47 tiles), trees (74 tiles), and water (1216 tiles).Figure 3
Large VHR dataset used and its qualitative results. (a) Original image of size 6500 × 10400 with example training tiles (10 per class), (b) testing map (ground truth), (c) legend, (d) NN classification result (OA = 89%), (e) final map (OA = 92.2%), and (f) tiles corrected by semisupervised learning.
(a)
(b)
(c)
(d)
(e)
(f)The second image is of size 14000 × 14000 and is taken from the city of Riyadh by GeoEye sensor in June 2010 (see Figure4(a)). The image has three spectral bands with a spatial resolution of 0.5 m. It has been divided into 4900 tiles of which 3860 are unlabelled, while 1040 are labelled into five land-cover types including bare soil (163 tiles), freeway (185 tiles), parking (109 tiles), residential (486 tiles), and sparse residential (97 tiles).Figure 4
Large VHR dataset used and its qualitative results. (a) Original image of size 13440 × 13440 with example training tiles (10 per class), (b) testing map (ground truth), (c) legend, (d) NN classification result (OA = 91%), (e) final map (OA = 96.5%), and (f) tiles corrected by semisupervised learning.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.2. Experimental Setup

Both large images have been subdivided into a grid of equal size tiles. The Jeddah image is divided in tiles of size 128 × 128 pixels, while the Riyadh image tiles are of size 196 × 196 pixels. We start the division from the top left corner, and any leftover pixels (in case the size of the large image is not a multiple of 128) appear at the right and bottom sides of the image. The ground truth classification maps are created manually by inspecting a random set of tiles (the test set) one by one and assigning one of the 9 labels to all tiles in the test set. Testing maps for both the Jeddah and Riyadh images are shown in Figures5(b) and 6(b), respectively. For both images, we extract feature representations from all tiles including HOG, LBP, BOW, and several pretrained CNN architectures. The set of pretrained CNNs used includes the imagenet-vgg-very-deep-16 and imagenet-vgg-m models by the VGG group at the University of Oxford [33, 34], the GoogleNet model developed by a team from Google Inc. [35], and Microsoft’s ResNet-152 model [36]. These CNNs are selected because they have won the first and second places in the ImageNet Large-Scale Visual Recognition Challenge- (ILSVRC-) 2014 and 2015 challenges. Table 1 shows the sizes of each feature representation.Figure 5
Results of output regularization via spatial graph Laplacian with sensitivity analysis of the parameterλ2 in (11), (a) for Riyadh dataset and (b) for the Jeddah dataset.
(a)
(b)Figure 6
Results of embedding a graph Laplacian in the hidden layer with sensitivity analysis of the parameterλ1 in (10). (a) Riyadh dataset with spectral graph. (b) Jeddah dataset with spectral graph. (c) Riyadh dataset with spatial graph. (d) Jeddah dataset with spatial graph.
(a)
(b)
(c)
(d)Table 1
Parameters used by various feature extraction techniques.
Feature type
Feature vector length
Note
HOG
4365
LBP
256
8 neighbourhood, uniform pattern
BOW
4000
Vocabulary size = 4000
CNN—vgg-very-deep-16
4096
CNN—GoogLeNet
1024
CNN—vgg-m
1024
CNN—ResNet-152
2048As for the assessment method, we use the Overall Accuracy (OA) and Average Accuracy (AA) measures. To assess the statistical significance of the improvements achieved by applying the semisupervised learning, we computed McNemar’s test, which is based on the standardized normal test statistic:(14)Zij=fij−fjifij−fji,where Zij measures the pairwise statistical significance of the difference between the accuracies of the ith and jth classifiers. The variable fij stands for the number of samples classified correctly and wrongly by the ith and jth classifier, respectively. Thus, fij and fji are the counts of the classified samples on which the two classifiers disagree. At the commonly used 5% level of significance, the improvement of classifier j over i is said to be significant if Zij<−1.96. Finally, all experiments are repeated 10 times with different random training tiles and then the average accuracy values are presented.


## 3.3. Comparison of Various Pretrained CNN and Finding the Best Hidden Layer Size

In this first experiment, we compare the performance of feature representations obtained by the pretrained CNNs and the ones obtained by BOW, LBP, and HOG. We also investigate in this experiment the sensitivity of the results with respect to the hidden layer size for all the types of feature descriptors considered in this study. We train a fully connected NN with the following hidden layer sizes 32, 64, 128, 256, 512, and 1024. We perform the experiment 10 times; each time, we select a different set of training sample tiles per class randomly and average the accuracy results.Figure7 shows the overall accuracies (OA) achieved for each type of feature descriptor and with different hidden layer sizes. It is clear from this experiment that CNN features are quite superior to all other types of handcrafted features regardless of the hidden layer size. This confirms other results about CNN that is reported in the object recognition literature.Figure 7
Comparison between various feature descriptors in terms of overall accuracy (OA). (a) OA for the Jeddah-IKONOS2 dataset and (b) OA for the Riyadh-GEOEYE dataset.
(a)
(b)Furthermore, if we consider CNN type features only, VGG very-deep-16 pretrained CNN, in particular, provided the best OA. We can also conclude that in general the range [64, 256] constitutes the best values for the hidden layer size. For the remainder of this paper, the hidden layer size of 256 and the VGG very-deep-16 pretrained CNN will be used in all experiments.


## 3.4. Semisupervised Learning Results

In this set of experiments, we present the results of applying semisupervised learning via graph Laplacian. We first implement embedding a spatial graph Laplacian as an output regularizer based on the RW algorithm without embedding any graph in the hidden layer. This step has a regularization parameterλ2 as shown in (11). In Figure 5, we show the Z score achieved for different values of the regularization parameter in the range of 1…20. From these results, we learn that output regularization does not provide significant improvement for the Riyadh dataset as the Z score is greater than −1.96.However, the improvement for the Jeddah dataset is significant (<−1.96). We also learn that the best value for the parameterλ2 is around 4 to 6; thus, we decided to fix this parameter to 5 for optimal performance of the method.In the second experiment, we present the results of embedding a graph Laplacian in the hidden layer. We investigated embedding both a spectral graph Laplacian and spatial graph Laplacian. The cost function to be optimized, shown in (10), has a regularization parameterλ1. Figure 6 shows the Z score achieved for different values of the regularization parameter in the range of 0.04…0.15. From these results, we learn that embedding a spatial graph Laplacian in the hidden layer degrades the classification accuracy significantly as the Z scores in Figures 6(c) and 6(d) are all above −1.96. On the other hand, embedding a spectral graph Laplacian in the hidden layer does provide significant improvement for both datasets as the Z score is less than −1.96. As for the best value for the parameterλ1, it is clear, from Figures 6(a) and 6(b), that it is around the values 0.08 to 0.1; thus, we decided to fix this parameter to 0.1 for optimal performance of the method.Next, we present the results of the complete system shown in Figure1. The results for the Jeddah and Riyadh image are shown in Tables 2 and 3, respectively. These tables show the OA, AA, and per class accuracy for 5 different scenarios: (1) using the fully connected NN only, (2) using the NN followed by spatial regularization using RW algorithm, (3) using the NN with spatial graph Laplacian in the hidden layer, (4) using the NN with spectral graph Laplacian in the hidden layer, and finally, (5) using the NN with both spectral graph Laplacian in the hidden layer followed by spatial regularization using RW algorithm, which is the complete solution proposed in this paper. These tables also show execution times in seconds, which include both training and testing times.Table 2
Classification results for the Jeddah dataset using one hidden layer of size 256.
OA
AA
Time (s)
Z score
Class 1
Class 2
Class 3
Class 4
Class 5
Class 6
Class 7
Class 8
Class 9
VGG-m
NN
86.05
72.84
0.23
0.00
60.32
82.35
54.40
82.38
56.58
73.77
77.03
71.72
96.97
NN + RW
86.42
73.24
0.04
−1.23
60.32
81.76
53.60
84.45
53.68
76.56
77.57
73.91
97.28
Spatial graph in the hidden layer
58.24
36.42
19.86
20.96
17.30
61.57
18.40
8.36
21.86
52.29
36.49
34.22
77.33
Spectral graph in the hidden layer
87.75
77.41
25.29
−2.71
80.00
84.90
57.40
81.77
80.57
59.34
86.22
69.69
96.81
NN + spectral graph + RW
88.33
78.08
0.04
−3.59
80.64
84.71
56.80
83.05
81.25
61.97
85.95
71.41
97.01
GoogleNet
NN
84.48
73.23
0.23
0.00
71.75
87.45
55.20
71.89
75.86
44.81
83.78
72.35
96.00
NN + RW
85.21
74.39
0.04
−2.44
72.86
87.45
55.20
75.80
76.94
46.12
84.59
74.38
96.14
Spatial graph in the hidden layer
77.37
58.45
20.05
6.59
47.30
87.65
11.60
78.72
70.21
35.57
58.92
45.31
90.81
Spectral graph in the hidden layer
85.97
76.82
26.38
−2.57
73.65
91.57
60.60
81.10
87.82
39.07
91.08
70.94
95.55
NN + spectral graph + RW
86.46
77.76
0.04
−3.32
74.92
91.77
60.80
84.27
89.33
37.60
91.89
73.59
95.67
ResNet
NN
87.68
79.33
0.37
0.00
81.59
83.14
63.80
91.04
75.13
70.00
81.35
73.59
94.36
NN + RW
88.33
79.69
0.05
−2.59
81.27
84.51
62.60
91.28
75.49
70.22
82.70
73.91
95.26
Spatial graph in the hidden layer
78.82
62.59
33.94
10.27
46.03
93.34
33.80
82.38
61.40
64.86
39.73
52.97
88.78
Spectral graph in the hidden layer
88.65
82.96
38.38
−2.00
84.76
91.96
78.00
91.95
79.22
69.73
86.76
70.16
94.12
NN + spectral graph + RW
89.01
82.92
0.06
−2.63
84.29
92.16
76.60
92.44
79.64
69.67
87.30
69.53
94.70
Very-deep-16
NN
88.41
78.52
0.62
0.00
74.76
88.63
67.60
96.34
57.67
76.45
82.16
65.63
97.45
NN + RW
89.37
78.98
0.09
−3.77
80.95
88.24
84.00
93.90
50.26
86.89
81.08
51.56
96.85
Spatial graph in the hidden layer
68.79
56.74
60.98
16.15
44.60
67.06
29.60
54.15
67.41
48.64
72.16
49.06
77.97
Spectral graph in the hidden layer
90.28
83.45
96.01
−3.00
82.86
90.59
74.80
93.96
74.77
73.28
95.67
68.28
96.86
NN + spectral graph + RW
91.23
83.68
0.09
−4.67
84.13
90.59
72.40
94.39
74.51
74.48
95.13
69.22
98.25Table 3
Classification results for the Riyadh dataset using one hidden layer of size 256.
OA
AA
Time (s)
Z score
Class 1
Class 2
Class 3
Class 4
Class 5
VGG-m
NN
87.96
81.98
0.22
—
99.09
71.49
86.26
96.53
56.55
NN + RW
88.07
82.21
0.05
−0.49
99.15
71.37
87.27
96.56
56.67
Spatial graph in the hidden layer
81.52
71.48
24.00
6.25
81.90
74.57
78.79
94.33
27.82
Spectral graph in the hidden layer
93.12
88.98
30.36
−5.30
99.22
85.03
91.01
98.61
71.04
NN + spectral graph + RW
93.18
89.07
0.05
−5.38
99.22
85.09
91.52
98.64
70.92
GoogleNet
NN
89.03
88.63
0.37
—
94.51
85.89
82.22
89.48
91.03
NN + RW
89.43
89.05
0.08
−1.06
94.51
85.43
83.64
90.09
91.61
Spatial graph in the hidden layer
80.97
75.54
23.97
5.76
96.21
77.77
70.61
85.50
47.59
Spectral graph in the hidden layer
93.45
91.03
55.17
−4.80
97.91
90.34
87.98
96.26
82.64
NN + spectral graph + RW
93.84
91.40
0.09
−5.20
97.91
90.40
88.18
96.83
83.68
ResNet
NN
89.58
86.33
0.35
—
88.30
81.09
87.78
95.40
79.08
NN + RW
89.72
86.44
0.07
−0.50
89.22
81.26
86.97
95.46
79.31
Spatial graph in the hidden layer
73.64
65.76
39.26
8.61
99.68
30.91
90.20
87.19
20.81
Spectral graph in the hidden layer
95.13
93.56
45.94
−5.86
98.76
89.54
96.36
97.50
85.63
NN + spectral graph + RW
95.29
93.80
0.07
−6.00
99.09
89.43
96.47
97.61
86.44
Very-deep-16
NN
90.97
88.09
0.63
—
95.36
75.60
94.85
96.81
77.82
NN + RW
93.13
91.84
0.11
−1.71
97.39
82.29
95.96
96.22
87.36
Spatial graph in the hidden layer
79.73
67.40
69.42
6.74
80.65
62.29
83.64
97.19
13.22
Spectral graph in the hidden layer
94.72
92.21
91.27
−4.72
98.00
88.85
95.12
98.29
80.79
NN + spectral graph + RW
95.26
92.91
0.11
−5.40
98.50
90.06
94.45
98.55
82.99The results show clearly that semisupervised learning has improved the accuracy significantly for all types of CNN. However, in agreement with an earlier experiment, the best results achieved are with thevgg-very-deep-16 pretrained CNN, with an OA of 91.23% and an AA of 83.68%. Also, the semisupervised learning has significantly improved the accuracy of the NN as the Z score reached −4.67. Similarly for the Riyadh image, the results shown in Table 3 indicate an even better improvement induced by semisupervised learning with a Z score below −5 for all types of CNN. Again, the best results are achieved using the vgg-very-deep-16 pretrained CNN, with an OA of 95.26%, an AA of 92.91%, and a Z score of −5.40.Finally, qualitative results are also shown in Figures3 and 4 for the Jeddah and Riyadh datasets, respectively. Here, we show the classification maps that have achieved the highest OA after semisupervised learning among the 10 successive runs of the algorithm. In Figures 3(e) and 4(e), we show the tiles whose labels are corrected by semisupervised learning. Comparing Figure 3(d) to Figure 3(c), one can see the type of improvement produced by semisupervised learning, namely, the reduced noise particularly in the water and in residential classes. As for the Riyadh classification maps shown in Figures 4(c) and 4(d), the improvements are clearly visible in the freeway class, where it is more recognizable as roads.In this final experiment, we compare the proposed method to a state-of-the-art method in the literature called the Laplacian extreme learning machine (LapELM) method [37]. To make the comparison fairer, we have changed the features used in the original method [37] to the same CNN features used here. Thus, the two methods now use the same CNN features extracted with the help of the vgg-very-deep-16 pretrained CNN. As for the other parameters of the LapELM method, we have set the ranges for the regularization parameter C to [0.001, 1000] and the RBF kernel parameter γ to [0.01, 5]. To estimate the best parameter values in these ranges, we use the deferential evolution technique, which outperforms the grid search technique as shown in [37]. As for the proposed semisupervised method, we note that the batch size parameter is usually set to 100 in the literature. However, here, we should take care not to set the batch size to be larger than the total number of training samples (number of samples per class × the number of classes). For this experiment, we set it equal to the number of samples per class times the number of classes.We also study the sensitivity of the method to reduced training samples (less than 10). Figure8 shows the OA for sample sizes from 2 to 10 for both methods being compared. Figure 8(a) shows the results for the Jeddah dataset, while Figure 8(b) shows the result for the Riyadh dataset. The first observation is that the classification accuracy for both methods degrades gracefully when the number of training samples per class decreases. However, the results confirm again the superiority of the proposed method compared to Laplacian ELM.Figure 8
Comparison between the proposed semisupervised NN and the Laplacian ELM method. (a) Results for Jeddah dataset. (b) Results for Riyadh dataset.
(a)
(b)


## 4. Conclusions and Future Work

In this letter, we proposed a practical and efficient solution to the semisupervised classification of large-scale VHR RS images. The proposed method is based on subdividing the image into square tiles and then extracting feature descriptors on these tiles. These in turn are fed into another fully connected neural network with semisupervised learning, which exploits the structural information contained in the unlabelled tiles. Semisupervised learning is implemented in our proposed method in two ways, first by embedding a graph Laplacian in the hidden layer and second by spatial regularization of the output layer using a random walker algorithm.Experimental results have shown that our proposed method can provide coarse classification maps for very large RS images with high accuracy. Furthermore, we have shown that our proposed semisupervised learning that combines a graph Laplacian in the hidden layer and spatial regularization of the output layer always provides significant improvement in terms of the overall accuracy, while having reasonable computational times for such large VHR images.As a future research direction, we can employ a self-paced selection strategy of unlabelled titles similar to the techniques used in [38, 39] to further improve the performance. One can also study different fusion techniques of the feature descriptors of many pretrained CNNs. Another direction is to research ways to combine data from different datasets in order to improve the accuracy of the learning algorithm.



---

*Source: 6257810-2018-04-05.xml*