2.1.

Deep Convolutional Neural Networks

Recent research has shown that neural networks trained on general image categorization tasks can be used as feature extractor for aerial imagery,^8,17 even though the viewing perspectives of overhead imagery and general images are quite different. This has allowed scene classification,¹² vehicle detection,¹¹ tree species mapping,¹⁸ and road detection¹⁹ in aerial images. The algorithms are applied on patches or tiles of the large aerial images as opposed to the entire image.

In deep convolutional neural networks (DCNNs), downsampling operation, e.g., max pooling, is performed after convolutional layers to capture the longer range contextual information and extract more abstract features. The downsampling process usually leads to coarse labeling or prediction at a lower spatial resolution. Fully convolutional neural networks (FCNNs) have been proposed to improve the coarse classification results and achieve dense end-to-end prediction. In Ref. 20, the authors proposed a skip architecture that combines semantic information from a deep, coarse layer with appearance information from a shallow, fine layer to produce accurate and dense semantic segmentation. Figure 3 shows the detailed architecture of FCN-8s skip network. In Ref. 21, authors proposed a semantic pixel-wise segmentation method called SegNet that utilizes deconvolutional decoder layers to map low-resolution encoder features to full input resolution features. The authors in Ref. 22 utilized an atrous convolution method to expand the support of the filter and removed most of the down-sampling operations to obtain dense labeling. These approaches have been successfully adapted to the dense semantic labeling of remote sensing images^5,6,23 and outperform the traditional pixel-level classification methods (such as SVMs²⁴) that employ hand-crafted feature descriptors. In our paper, we adopt the FCN in Long et al.²⁰ paper as our pretrained network and fine-tuned it using the VHR optical imagery with a two-stage training scheme. We utilized the per pixel probability output from the network as one of the local potential functions in our higher-order CRF framework.

Fig. 3

Proposed FCN-8s architecture with two skip layers. This network architecture learns to combine coarse, deeper layer information with fine, shallow layer information. It combines predictions from the final layer and the pool4 layer to generate a finer prediction. Including the output from pool3 can provide an even more detailed prediction.

2.2.

Conditional Random Fields

In CRF models, each pixel is assigned a label by examining its color or spectral characteristics and also spatial relations. Pixels of an object can have complex local or global dependencies. For image classification, in general, two important assumptions are made with regard to spatial smoothness and contextual coherence: (1) neighboring pixels tend to belong to the same class except on the object boundaries. (2) adjacent pixels/objects from natural images must follow a certain practical meaning, e.g., the car is more likely found on the ground than on a tree. CRFs can be used to model these dependencies at a local and global scale. These energy-based random field methods have been widely used for exploiting contextual information for both computer vision and remote sensing images.^{1–3,25–28} The full potential of the random field methods has not been realized due to the limitation of inference methods for higher-order node connections, particularly for the large-scale data. Recently, research has been done to employ higher-order random fields with efficient inference methods.^29–33 The robust higher-order potential and graph cut inference method^29,30 has stood out due to its efficiency on the relatively large-scale dataset and the state-of-the-art semantic segmentation performance. It has been applied to the remote sensing applications such as road and rooftop extraction.^25,26

2.3.

Multisensor Classification

There are two commonly used fusion techniques for the classification of multisensor images. In one methodology, features are initially extracted from all imaging modalities and then fused together by concatenation or selection. The combined features are then utilized in a supervised training scheme to obtain a label map.³⁴ In the second fusion method, classification results are initially obtained for each imaging modality, and then these predictions are merged to achieve an optimal output (referred to as decision-level fusion^35–37). Sherrah⁵ trained two separate neural networks for each modality and concatenated the learned features at the last convolutional layer. In Ref. 7, multiple predictions from individual modalities are averaged. They also introduce a network that combines features at the final layers by a residual correction framework. The later network achieved slightly better quantitative results in comparison with the averaging of outcomes. Recent research in Audebert et al.³⁸ provides a solid investigation of applying early and late fusion in deep neural networks for LiDAR and multisepctural data.

Authors in Ref. 39 represent all the bands of multisensor data as a single image and utilize it to train a multiresolution CNN. They also extract hand-crafted features from the image to train a second classifier. They combine the individual results at a decision-level. Our proposed method differs from the previous algorithm in three different aspects. First, we apply fully convolutional neural networks only on the three bands (optical imagery) of multisensor data, whereas in their work, the input is a combination of all bands (requires training of neural network from random initialization). As our network uses only three bands, pretrained network weights could be fine-tuned with existing limited ground truth. Next, we learn the weight parameters to combine individual probabilistic outputs in a CRF framework using training data. Their approach combines results from two classifiers directly in a rigid way. Lastly, our method uses a higher-order CRF model, whereas their model uses only standard unary and pairwise potentials.

3.1.

Deep Fully-Convolutional Neural Network

VHR optical imagery provides rich low-level and high-level features. To fully take advantage of such information, we used features from neural networks to generate high-resolution per class probability prediction. Pixel-wise semantic segmentation of an image through neural networks was initially proposed by Long et al.²⁰ by extending a popular image classification neural network named VGG-16. The VGG-16 network consists of multiple convolutional layers followed by two fully connected and a final scoring layer. The last layer acts as a classifier to give scores to different predefined objects. The authors in Ref. 20, in their work, modified the fully connected layers to add spatial support and thereby generate the labels for each pixel. Specifically, the output of the fully connected layer has dimensions of $n\times n\times 4096$ instead of $1\times 4096$, where $n$ indicates the size of the down sampled spatial support of the fully connected layer. They called the resultant network the fully convolutional neural network (FCN). FCN-8s contains five layers with multiple convolutions and rectified linear activation functions. Each of these layers is succeeded by a max pooling (downsampling) operation. The architecture is shown in Fig. 3. The sixth and seventh layers are two fully connected convolutional layers. The eighth convolution (score) layer generates outputs corresponding to the number of classes in the ground truth. Upsampled outputs from the eighth layer are combined with the outputs from pool 4 layer. The result is again upsampled and merged with pool 3 layer outputs. Fusing the output of pool 3 and pool 4 layers (skip connection) assists in obtaining finer semantic labels. The fine-tuned network was then used to generate per pixel probability maps, which are denoted as $P_{\mathrm{FCN}}$.

One of the FCN-8s classification results and its corresponding error map are shown in Figs. 4(c) and 4(g). We notice that part of the rooftops is misclassified as impervious surface. This demonstrates that objects and land covers sometimes can have very similar spectral and spatial characterizations in aerial imagery and thus FCN-8s is not able to correctly distinguish them. Fortunately, LiDAR data provide us the complementary information, i.e., elevation, to improve the performance of classification.

Fig. 4

Classification results of ISPRS Potsdam dataset 4_15 and their corresponding error maps. (a) Color-infrared images, (b) classification results of MLR using LiDAR features, (c) classification results of FCN-8s using IR, R, and G channels, (d) classification results of our proposed method, (e) normalized DSM map, (f) segmentation map obtained from the segmentation algorithm. (g) and (h) Corresponding error maps for (c) and (d), respectively, where red pixels indicate misclassifications.

3.2.

Multinomial Logistic Regression

The LiDAR data are provided as normalized digital surface maps (nDSMs), which contain limited contextual information in comparison with the VHR optical imagery. Previous works demonstrated that training an additional neural network on the LiDAR data can increase the overall performance of the classification task. We, however, hypothesize that utilizing an efficient linear classifier shall be sufficient to take advantage of the LiDAR information. In our work, we employ an MLR on hand-crafted features derived from both LiDAR data and optical imagery. MLR has been used in numerous remote sensing classification tasks,⁴⁰ and it can produce multiclass probabilistic estimates. More importantly, MLR is fully probabilistic, therefore it provides calibrated probabilities off-the-shelf, whereas SVM and random-forest require postprocessing to compute multiclass (e.g., one versus all) and probabilities (Platt method⁴¹ using cross validation). The probability that a pixel $x_i$ belongs to a particular category $c$ can be calculated as follows:

3.3.

Higher-Order Conditional Random Field Formulation

The two probabilistic predictions generated using the optical and LiDAR data, respectively, provide complementary information to the optimal class labeling. FCN-8s trained on the visible bands exploits more spectral and spatial information and, therefore, it learns object structure much better than MLR does. For land covers/objects that lack textures or have indistinguishable spectral characterization, MLR takes advantage of height information and generates very reliable and efficient predictions. We combine these two outputs in the higher-order CRF framework such that the final classification results can have spatial/contextual coherence and in addition preserve object boundaries. The CRF is an undirected graphical model defined as follows:

3.3.3.

Higher-order potential

The higher-order potential extends the smoothness assumption from the neighboring pixels to the local regions. The idea is to first group pixels in the images into coherent regions so that each region ideally belongs to one object. Using an edge awareness segmentation algorithm, this higher-order potential is particularly useful in preserving the object boundaries. Given a set of segments denoted as $C$, we can define our higher-order potential as

Eq. (13)

$${\psi }_h(x)=\sum _{c\in C}{\psi }_c(x_c).$$

Segments $C$ are usually generated by an unsupervised segmentation algorithm.^13,45 We choose to take the form of the robust $P^N$ Potts’ potential proposed in Ref. 29 for ${\psi }_c(x_c)$, which has been proved to be particularly useful. It takes the form of

Eq. (14)

$${\psi }_c(x_c)=\{\begin{array}{ll}{N}_i(x_c)\frac{1}{Q}{\gamma }_c^{\max },& \text{if}{N}_i(x_c)<Q\\ {\gamma }_c^{\max },& \text{otherwise}\end{array},$$

where $N_i(x_c)$ denotes the number of pixels taking different labels from the dominant label of the segment. ${\gamma }_c^{\max }$ is the maximum cost that will be added. ${\gamma }_c^{\max }$ can be defined in a way that takes into account the quality of the segmentation and the contextual information. We define ${\gamma }_c^{\max }$ in our case as the following:

Eq. (15)

$${\gamma }_c^{\max }={\theta }_c|c|·e^{[-{\theta }_h\frac{\sum _{i\in c}{(I_i-{\mu }_i)}^2}{|c|}]},$$

where $|c|$ counts the number of pixels in the segment $c$. $I_i$ is the features of each pixel in the segment. $\mu$ is the mean observations over the current segment. The exponential term in Eq. (15) calculates the variance of observations within one segment, which indicates the quality of the segment. Therefore, the higher-order potential will impose the cost accordingly. For instance, if the variance is small, the segment very likely captures a coherent region. Therefore, taking different labels in those coherent regions can cost more penalties. On the contrary, when the variance is large, the cost of taking different labels in one segment is small. That is because large variance usually means that the region might contain multiple objects so that forcing all the pixels in this segment to take the same label can cause misclassification. Unlike the $P^N$ Potts model, which imposes a strict cost to any heterogeneous labeling within one segment, the robust $P^N$ Potts model proposes a linear truncated cost that is dependent on the number of pixels taking different labels from the dominant label of the segment. The heterogeneity of the labeling is controlled by the parameters ${\theta }_c$ and $Q$, respectively.

Now let us combine the Eqs. (5), (6), (10), and (13) to formulate our higher-order CRF energy function as

Eq. (16)

$$E(x)=\sum _{i\in v}{\psi }_i(x_i)+\sum _{i\in v,j\in N(i)}{\psi }_{ij}(x_i,x_j)+\sum _{c\in C}{\psi }_c(x_c),$$

where ${\psi }_i(x_i)$, ${\psi }_{ij}(x_i,x_j)$, and ${\psi }_c(x_c)$ take the form in Eqs. (7), (12), and (14), respectively. Figure 5 shows the configuration of our proposed higher-order CRF. To be noticed, for the global node of a segment, it can be one of the predefined labels in $L^M$, or if there is no dominant label for the segment, the segment can take a free label $L^F$, which means that CRF imposes no cost on heterogeneous labels in this segment. In such a case, it is equivalent to only applying the pairwise CRF.

Fig. 5

Probabilistic graphical model for the proposed higher-order CRF.

3.4.

Learning Conditional Random Field Parameters

The parameters of the proposed higher-order CRF are learned in a stepwise training procedure. We first learn to combine the two probabilities obtained from the FCN-8s: $P_{\mathrm{FCN}}$ and MLR: $P_{\mathrm{MLR}}$ to form the unary potential ${\psi }_i(x_i)$. We construct the unary potential using a softmax classifier to recast the two probabilities into a single probability. The parameters ${\sigma }_m$ can be learned using a cross-validation set with a standard maximum likelihood estimation (MLE) procedure. We then keep ${\sigma }_m$ constant and proceed to train the pairwise CRF parameters: ${\theta }_{\alpha }$, ${\theta }_{\beta }$, ${\theta }_{\gamma }$, and the matrix $T(x_i,x_j)$ without the higher-order term using an approximate marginal inference method proposed in Ref. 46. The marginal inference procedure takes into account model mis-specification and inference approximation. Finally, the higher-order parameters such as ${\theta }_c$ and $Q$ are learned by performing cross validation within an empirical range. Moreover, the value of ${\theta }_h$ is picked empirically, and we use ${\theta }_h=2$ for all of our experiments. The tuning of parameters ${\theta }_c$ and $Q$ will be discussed in Sec. 4.

3.5.

Inference Using Graph Cuts

The inference problem in a higher-order CRF is to find the set of class labels that minimizes the proposed energy function in Eq. (16), i.e., $\arg {\min }_{x}E(x)$. In general, this energy optimization problem is an NP-hard problem. However, there are certain classes of tractable functions that can be solved in polynomial time, e.g., submodular functions. The energy that we propose happens to belong to those functions. There are two main methods to approximate such an energy minimization problem: message passing algorithms such as belief propagation, and “move making” algorithms, such as graph cut. In this paper, we choose to use the “move making” graph cut algorithm due to its computational efficiency.

The “move making” graph cuts inference algorithm (e.g., $\alpha$-expansion and $\alpha \beta$-swap) has been successfully used to infer the higher-order CRFs.^29,30 We will review the $\alpha$-expansion graph cut and then show how to deal with the higher-order CRFs with $\alpha$-expansion graph cut by adding auxiliary nodes. We refer readers to papers^29–31 for a more detailed description.

The “move making” graph cut algorithm was proposed to efficiently solve the multiclass classification using s–t min-cut-based graph cut algorithm, which was first introduced for binary classification problems.⁴⁷ As shown in Fig. 6(a), the “move making” algorithm usually starts from an initial set of labels and then iteratively updates the labels to find the solution that has the lowest energy. To utilize s–t min-cut algorithm, each update in the “move making” algorithm has to be a binary decision. For instance, $\alpha$-expansion allows any random variable to either maintain its current label or take the proposed label $\alpha$. We form a transformation function that converts the energy from the label space into the “move” space. We then deduce its corresponding move energy. The transformation function $T_{\alpha }(·)$ for the $\alpha$-expansion as

Fig. 6

(a) Flowchart of the move-making graph cut algorithm, (b) S–t min-cut for $\alpha$-expansion, and (c) $\alpha$-expansion for higher-order CRF by adding auxiliary variables.

3.6.

Segmentation for Higher-Order Potential

In the higher-order CRF framework, the quality of the segments impacts the success of the final labeling. It is critical to choose a robust segmentation algorithm that can yield a dense semantic labeling with fine boundaries. GSEG algorithm is one such unsupervised multichannel image segmentation algorithm that utilizes the gradient histogram acquired from the color images to iteratively cluster pixels from lower gradient to higher gradient. The GSEG algorithm is primarily based on color-edge detection, dynamic region growth, and a unique multiresolution region merging procedure. The detailed description of the algorithm can be found in Ref. 13.

As GSEG uses the gradient histogram, which helps to preserve the object boundaries. See one of the illustrations of GSEG segmentation results with different initial seed sizes in Fig. 7. GSEG is particularly suited for aerial image segmentation because the size of the objects in the VHR aerial images can vary from a several hundred pixels to tens of thousands. GSEG does not pose strict constraints of the segment size. Instead, GSEG keeps growing coherent regions until no local low-gradient pixels can be found. It is, therefore, able to generate segments for objects different scales. The performance of GSEG is mainly impacted by two parameters, one is the initial seed size $\tau$, and the other is similarity ratio $\gamma$. The former is to determine the initial size of the low-gradient region, and the latter controls the sensitivity of local feature change.

Fig. 7

(a) RGB image patch from Potsdam dataset. (b), (c), and (d) Segmentation results of GSEG algorithm with initial seed size of 15, 100, and 150, respectively.

4.1.

Remote Sensing Datasets

Potsdam has a ground sampling distance of 5 cm and includes optical orthophotos with four spectral channels, IR, R, G, B and the corresponding coregistered normalized DSMs. The entire collection is divided into 38 image patches, from which 24 images with ground truth labels are used for training, and the remaining 14 images are employed for testing. The Vaihingen dataset has a ground sampling distance of 9 cm and consists of 33 images, with 17 images for training and 16 for testing. The optical images have three spectral channels only: IR, R, and G are accompanied by the coregistered DSMs. The imaging data of Zeebruges were acquired using an airborne platform flying at the altitude of 300 m over the urban and harbor areas of Zeebruges, Belgium (51.33 N, 3.20 E). The data were simultaneously collected and georeferenced to WGS-84. The point density for the LiDAR sensor was $\sim 65\text{points}/{\mathrm{m}}^2$, which translates to a point spacing of $\sim 10\mathrm{cm}$. The color orthophotos were taken at nadir and had a spatial resolution of about 5 cm. The dataset is organized into seven separate tiles. Each tile includes a $\mathrm{10,000}\times \mathrm{10,000}\text{pixel}$-sized portion of the color orthophoto (GeoTIFF, RGB), and a corresponding $5000\times 5000$ DSM. Five of the tiles are used for training, and the remaining two for testing.

In our work, we selected the following five classes: impervious surface, buildings, low vegetation, trees, and cars from the Potsdam and Vaihingin given that their ground truth is readily available. For Zeebruges, we selected two additional categories of interest: boats and water for the same reasoning. Figure 8 shows one of the training images in Zeebruges along with its DSM and ground truth. Figure 9 shows the color legend for each class and the corresponding categories that were identified in each dataset.

Fig. 8

Example of training image from the Zeebruges dataset. (a) RGB orthophoto, (b) digital surface model, and (c) ground-truth labels.

Fig. 9

Semantic labeling color legend and the corresponding categories that are classified in each dataset.

4.2.

Training the Fully-Convolutional Neural-8s Network

To effectively train and test our proposed algorithm, we selected an image resolution of $224\times 224$ for both ground truth and test images, respectively. The image patches could be cropped from the dataset in numerous ways. In our experiments, we generated the training set using the following guidelines: (a) for each class, randomly select nonoverlapping image patches of size $224\times 224$ from each image, (b) for each class, randomly select 1000 pixels in each training image and obtain a $224\times 224$ patch with selected pixel as starting point, and (c) randomly choose 50 cars from the training images to ensure adequate car samples (and also select additional boats in Zeebruges dataset). The process outlined above yielded a training image set of 43,516, 18,780, and 36,000, using the Potsdam dataset, Vaihingen dataset, and Zeebruges dataset, respectively.

We initialized the parameters of the FCN-8s network with pretrained weights, which were obtained using a large dataset of color images and corresponding labels.²⁰ The weights were then fine-tuned using the corresponding training dataset.

The training process includes two stages. In stage 1, we initialized the weights for convolutional layers 1 through 5, and fully connected layers 6 and 7 with pretrained weights that were kept constant during the first learning phase. We then train the remaining layers—with the exception of the connection and score layers—using randomly initialized weights. This training was done with a learning rate of $1\times {10}^{-3}$ for 35 epochs. The learning rate was decreased by 0.1 after 15 and 30 epochs. In stage 2, we set the parameters using the learned weights from the previous stage and then fine-tuned all of the layers with a reduced learning rate of $1\times {10}^{-5}$ for 35 epochs. Again, the learning rate was decreased by 0.1 after 15 and 30 epochs. A stochastic gradient descent algorithm was utilized for the fine-tuning using the Caffe [Caffe: Convolutional Architecture for Fast Feature Embedding] toolbox.⁵¹ For the test images, and to avoid blocky artifacts, we chose tiles with a stride of 112, i.e., with an overlap rate of 50%. The tiles were then processed by the network for classification.

Our method only requires training one network compared with other competitive deep learning-based fusion frameworks, such as Refs. 5 and 52. The number of parameters for training one FCN-8s is 134,279,076 and for two FCN-8s is 268,558,152. The training time for one FCN-8s network takes 8 h and 15 min for the first 35 epochs and 20 h and 30 min for the second 35 epochs. While for training two FCN-8s, it takes 10 h 49 min for the first 35 epochs and 20 h and 45 min for the second 35 epochs. Training one neural network consumes much less memory as well as takes less training time.

Fig. 10

Semantic labeling results of our proposed higher-order CRF method with different segments. (a) RGB image patch from Potsdam training set, area 4_10, (b) ground-truth labeling, (c) results of FCN-8s trained on IR, R, and G channels only. (d) and (e) Results of our proposed higher-order CRF with GSEG algorithm using initial seed size of 15 and 100, respectively. (f) Classification result of our proposed higher-order CRF with ground-truth segments (upper bound).

4.3.

Learning the Higher-Order Conditional Random Field Model

The classification results of the higher-order CRF are, in general, affected by the quality of the segmentation map,^5,29 which is employed by the algorithm. To minimize misclassification, we utilized the GSEG algorithm¹³ to provide all underlying segmentations due to its optimal performance against the state of the art.⁵³ To this effect, we have generated multiple segmentations using various parameters to test the robustness of our proposed higher-order CRF algorithm. The results were compared against the available ground truth that was generated from the training dataset, where each segment articulates the proper object boundaries.

We found that different segmentations do affect the labeling performance of our higher-order CRF fusion method. At the initial seed size of 15 pixels, the vehicle $F_1$-score is even lower than the one without using CRF and fusion of LiDAR (see the comparisons in Table 1, where HCRF_# means that it uses our proposed higher-order CRF framework, and the segments are obtained by using GSEG algorithm with initial seed size of #). With an increase in the initial seed size, the vehicle $F_1$-score improves a noticeable amount, and every other category’s $F_1$-score keeps increasing and peaks at the initial seed size of 100. Comparing with using the ground truth segmentation, our proposed method achieves an overall accuracy of 96.05%. This demonstrates that our proposed higher-order CRF can leverage the global knowledge if given correctly. Although we acknowledge that ground-truth segmentation map is rarely accessible by an unsupervised segmentation technique, we would argue that by choosing a suitable segmentation algorithm with its appropriate parameters, applying higher-order CRF tends to improve the final dense semantic labeling (Fig. 10).

Table 1

Quantitative results of three validation images with different segments. FCN-8s_CIR: fully-convolutional neural network trained on only IR, R, and G channels; nDSM: multisensor fusion with normalized DSM; PCRF: pairwise CRF; HCRF_#: higher-order CRF using segments with initial seed size of #. HCRF_GT: ground truth segmentation, used as our upper-bound performance.

Table 2

Experimental comparison of using two different segmentation algorithms: SLIC45 and GSEG. In addition, we also compared the fusion of FCN-8s and the prediction with LiDAR features using SVM classifier and MLR.

4.4.

Results of the Potsdam Dataset

We tested our proposed method on the Potsdam dataset and compared its performance qualitatively and quantitatively against the most competitive methods that are published on the ISPRS benchmarking website. The results are shown in Fig. 11, and tabulated in Tables 1 and 3, respectively. SVL_1 serves as the baseline method and does not employ a neural network-based design.⁶ UZ_1 utilizes a CNN-based system relying on a downsample-then-upsample architecture to generate dense semantic labeling without any postprocessing. KLab_3 adopts a DCNN for classifying multisensor remote sensed images. DST_6 employs two separate FCNs for optical and LiDAR data, respectively, incorporates a pairwise CRF at the end of their FCN framework to gain a 0.6% boost in overall accuracy and omits the downsampling operation that has been used in most FCNs.⁵ The qualitative results of our proposed method can be found in Fig. 11 (denoted as DNN_HCRF in Table 3 and RIT_L7 on the benchmarking website). The results tend to preserve more object boundaries when in comparison with the FCN-8s-based technique (see Fig. 11 column #2). This fact is also shown in Table 3 across the various classes. Compared with Ref. 5, our method eliminates the expense of training the second neural network and utilizes the higher-order CRF framework to improve the overall accuracy by taking advantage of the context information found in most images. This yields more refined object boundaries (especially for buildings) as shown in Fig. 11 and confirmed in Tables 1–3 by the quantitative accuracy. Most of the errors observed in our experiments are the result of inaccuracies in the underlying segmentation maps.

Fig. 11

Results on the Potsdam dataset. Rows I, II, and III illustrate testing patches 3_14, 5_15, and 7_13, respectively. Column #2 shows the results from FCN8s trained only on IR, R, and G channels and column #4 provides the results of our proposed method. The corresponding classification error maps are shown in column #3 and #5, where the red regions indicate the misclassified pixels.

Table 3

Results on 14 test images of the Potsdam dataset. DNN_HCRF: results of using proposed higher-order CRF fusion method. SVL_1, UZ_1, KLab_3, and DST_6 are the published benchmarking methods on ISPRS Potsdam 2-D labeling contest website.

4.5.

Results of the Vaihingen Dataset

The Vaihingen dataset has a slightly lower ground spatial resolution (9 cm) compared with the Potsdam dataset (5 cm). Furthermore, since Vaihingen dataset does not provide the corresponding normalized DSMs, these were generated by manually filtering the ground points. The qualitative and quantitative results—along with the most up to date benchmark comparisons from the ISPRS website—are shown in Fig. 12 and Table 4, respectively. It should be noted that our proposed technique performed adequately well as compared with state of the art. Errors were primarily due to the underlying segmentations and manually generated DSMs that contain large flat regions with gradually changing elevation. Those regions can be misclassified as buildings (e.g., see Fig. 12, row #2, column #5). Note that the Vaihingen dataset contains more flat areas with gradually changing elevation than in Potsdam dataset. Both the FCN-8s and MLR are unable, in general, to correctly recognize the gradually elevated regions without an accurate nDSM.

Fig. 12

Results on the Vaihingen dataset. Column #1 and #4 include four Vaihingen dataset test patches, namely area_6 (column#1 and row#1), area_12 (column#4 and row#1), area_22 (column#1 and row#2), and area_20 (column#4 and row#2). Column#2 and column#5 are the labeling results of our proposed method. Column#3 and column#6 are the error maps where the red regions indicate the misclassified pixels.

Table 4

Results on 17 test images of the Vaihingen dataset. SVL_3, ADL_3, ONE_7, UZ_1, DLR_10, and UOA are the published benchmarking results on the ISPRS Vaihingen 2-D labeling contest website. Our proposed method is denoted as DNN_HCRF in the paper or RIT_L7 on the benchmarking website.

4.6.

Results of the Zeebruges Dataset

The Zeebruges dataset was first introduced in the 2015 IEEE GRSS data fusion contest. The most competitive experimental results of this contest have been published in Ref. 10. The best numerical result reported in Ref. 10 is achieved by completely retraining a CNN end-to-end using R, G, B, and height data. Several non-CNN-based techniques are listed as the baseline methods. As shown in Table 5, the results using a pretrained FCN-8s outperform the best results reported in the 2015 contest. Our proposed method further improves the overall accuracy by a large margin (2.35%) compared with only applying FCN-8s. Figure 13 shows the visual improvements in object boundaries by using our proposed higher-order CRF. This is attributed to the high-resolution (ground resolution of 5 cm) RGB images and the extremely high density of LiDAR data. VHR aerial imagery in the Zeebruges results in a high quality of segmentation, which subsequently boosts the performance of our proposed higher-order CRF framework.

Table 5

Results on the Zeebruges dataset. HSVDGr/SVM, blesaux, and RGBd+ trained AlexNet are the most recent benchmarking results presented in the review paper.10 FCN-8s only uses FCN-8s without CRFs. DNN_HCRF is the proposed method.

Fig. 13

Qualitative results on the Zeebruges dataset. The illustration demonstrates that our proposed higher-order CRF can produce more accurate object boundaries compared with FCN-8s as seen from the zoomed images.

5.1.

Fine-Tuning of the Pretrained Neural Network

In addition to the training procedure introduced in Sec. 4.2, we also trained the FCN-8s neural network using a one-stage (namely DNN_T1) instead of a two-stage training strategy (namely DNN_T2). To this effect, we initialized the parameters of convolutional layers 1 to 5 and fully connected layer 6 and 7 with pretrained weights in Ref. 20 and assigned random weights for the remaining layers. The parameters of convolutional layers 1 to 5 were held constant during the training, and the parameters of the other layers were fine-tuned with a learning rate at $1\times {10}^{-3}$, 35 epochs, multiplying learning rate at 15 and 30, up by 0.1. We evaluated these two training strategies on the Potsdam testing dataset and found that the two-stage training strategy outperformed the one-stage training strategy regarding overall classification accuracy by 2.3% as shown in Table 6.

Table 6

Classification results using only FCN8s with two different training strategies. The one-stage training method is denoted as DNN_T1, and the two-stage training method is denoted as DNN_T2.

As our experiments showed, different training strategies for learning fully convolutional neural networks parameters have a significant impact on the overall accuracy of semantic labeling. We think that the second step of fine-tuning all the layers in the two-stage training strategy is attributed to the performance improvement. Because the pretrained neural networks are usually trained based on general viewing perspective of objects, the shallow convolutional layers also need to be fine-tuned to adapt to the different appearance of aerial images. We found that changing the fine-tuning strategy significantly impacts the classification results.

5.2.

Parameter Learning of the Higher-Order Term

Based on Eq. (18), to achieve the optimal performance of the higher-order CRF, we need to train the hyper-parameter $\frac{1}{Q}$. $\frac{1}{Q}$ can be interpreted as a truncation term that determines the percentages of the number of pixels in one segment are allowed to take a different label from the dominant label for the segment. We refer the readers to this paper¹⁴ for more details on the hyperparameters training process.

5.3.

Higher-Order CRF Modeling

We have discussed the quantitative improvements of using the higher-order CRFs compared to only using the pairwise CRFs in the previously published paper.¹⁴ We have seen the same quantitative and qualitative improvements for the Vaihingen and Zeebruges dataset. Especially for the Zeebruges dataset, as it has the even higher spatial resolution RGB images, the higher-order CRF gains more advantages in terms of resolving the object boundaries (see Fig. 13).

The need for a higher-order CRF is discussed in Ref. 39, in which the authors argued that a higher-order CRF sometimes had an adverse impact on classification accuracy. We agree on the point that the performance of higher-order CRFs is sensitive to the quality of segments, which is scene dependent. As shown in Table 1, the vehicle $F_1$-score drops when higher order CRF takes a small initial seed size, which resulted in over-segmented objects. Based on our experiments, as long as we choose a proper segmentation algorithm and find its appropriate parameters, using higher-order CRF gains an overall quantitative and qualitative improvement for dense semantic labeling compared with only using pairwise CRF. Furthermore, incorporating a higher-order CRF provides potential opportunities for further improvement by utilizing object-level contextual information in a hierarchical random field, as proposed in Refs. 3031.–32.

The main contribution of our paper is mainly focusing on the qualitative and quantitative improvements by introducing the higher-order CRFs after we obtained the probabilistic results from the state-of-the-art CNNs. The individual results from each CNN architecture may differ. As we treated these two parts (CNNs and CRFs) separately, we expect improvements even if we use other CNN architectures. Higher-order CRFs leverage the boundary and context information that is not necessarily learned by CNNs. The performance of our approach combining with other CNN architectures needs more investigation, but one would expect incremental improvement in accuracy by carefully comparing and selecting other deep learning architectures.