Satellite Image Super-Resolution via Multi-Scale Residual Deep Neural Network

Abstract

Recently, the application of satellite remote sensing images is becoming increasingly popular, but the observed images from satellite sensors are frequently in low-resolution (LR). Thus, they cannot fully meet the requirements of object identification and analysis. To utilize the multi-scale characteristics of objects fully in remote sensing images, this paper presents a multi-scale residual neural network (MRNN). MRNN adopts the multi-scale nature of satellite images to reconstruct high-frequency information accurately for super-resolution (SR) satellite imagery. Different sizes of patches from LR satellite images are initially extracted to fit different scale of objects. Large-, middle-, and small-scale deep residual neural networks are designed to simulate differently sized receptive fields for acquiring relative global, contextual, and local information for prior representation. Then, a fusion network is used to refine different scales of information. MRNN fuses the complementary high-frequency information from differently scaled networks to reconstruct the desired high-resolution satellite object image, which is in line with human visual experience (“look in multi-scale to see better”). Experimental results on the SpaceNet satellite image and NWPU-RESISC45 databases show that the proposed approach outperformed several state-of-the-art SR algorithms in terms of objective and subjective image qualities.

1. Introduction

Remote sensing satellites, which observe objects on the ground from outer space, are widely used in various real applications, such as environmental monitoring, resource exploration, disaster warning, and military applications. The observed images from satellites generally have low-resolution (LR) due to the limitations of spaceborne imaging equipment (Charge-coupled Device (CCD) sensors) and communication bandwidth. In addition, satellite images are affected by atmospheric turbulence, transmission noise, motion blur, and undersampling optical sensors. The quality and resolution of images from remote sensing satellites cannot meet the requirements of real satellite image analysis. Super-resolution (SR) technology can overcome hardware limitations and improve the spatial resolution of images through software manner. The first SR algorithm [1] was developed to improve the resolution of Landsat remote sensing images by fusing multi-frame complementary information. In the last decades, SR has been successfully applied to enhance the resolution and quality of remote sensing satellite images. A well known example is SPOT-5, which reaches 2.5 m resolution through the SR of two 5 m images that are sampled from shifting a double CCD array by subpixel sampling interval [2,3]. Traditional SR image generation methods usually require multiple spatial/spectral/temporal low-resolution images of the same scene [4,5].

Existing image SR algorithms are divided into two categories, namely reconstruction- and learning-based algorithms [6]. Reconstruction-based algorithms fuse subpixel LR multi-frame information and reconstruct their latent high-resolution (HR) images. Previous satellite SR methods utilize reconstruction-based methods in solving the inverse problem of the degradation process. Reconstruction-based methods model the degradation process of imaging with mathematical formulas by using degradation factors, such as downsampling, optical blur, atmospheric disturbance, registration error, geometric deformation, and motion compensation [7,8,9]. Although reconstruction-based methods are simple and intuitive and can be flexibly combined with prior constraints, they rely on accurate subpixel precision estimation.

Inspired by the immense success of machine learning in object recognition and other tasks, learning-based SR methods have been highly valued and have become the mainstream direction of research. They aim to learn a mapping function between LR and HR image/patches through the prior information provided by a training dataset. Learning-based SR algorithms can obtain better subjective and objective reconstruction performance than reconstruction-based methods because external training databases provide considerable a priori information. In terms of the usage of prior training samples, learning-based SR algorithms can be divided into three categories, namely regression-, representation-, and deep-learning-based algorithms. Some representative regression-based [10,11,12] and representation-based SR algorithms [13,14,15] yield decent subjective and objective performance. These methods are efficient with flexible framework for using regularization terms.

Deep-learning-based approaches provide an end-to-end solution for learning complex mapping functions and are rapidly and successfully applied on SR tasks. A complex nonlinear mapping relationship between LR and HR patches is learned through convolutional neural networks (CNNs) [16,17], considering their excellent learning capability. Shi et al. [18] constructed a subpixel CNN, which provides a novel manner of directly and efficiently learning the mapping function from LR to HR images, which is further efficient. Kim et al. [19] stated that the construction of a deep network can effectively alleviate training difficulty with deep residual learning. Lai et al. [20] built a pyramid network for fusing multi-scale residuals in the feature domain. A generative adversarial network (GAN) [21,22], which comprises generator and discriminative networks, was used to generate fake details for simulating a good visual output. For satellite images, Luo et al. [23] replaced zero-padding with self-similarity to avoid the addition of unusable information and achieved good results. Wang et al. [24] proposed a multi-memory CNN for video SR to retain inter-frame temporal correlations.

The above-mentioned SR approaches mainly focus on the general nature images. As for satellite image, the object scale in the image is relatively different due to wide-range imaging, and it has important roles in vision tasks, such as segmentation, feature extraction, and object tracking. Some deep-learning algorithms designed for general images cannot efficiently handle satellite images because they do not specially consider the multi-scale nature of satellite images. Moreover, adequate high-frequency information, such as edges and textures, are crucial for satellite image detection [25] and object recognition [26,27,28,29,30]. The use of a single structure network in predicting and reconstructing objects without considering their different scales results in poor reconstruction performance. One practical solution is to explore the multi-scale information into deep neural networks. Zhang et al. [31] used multi-scale spatial structural self-similarity to learn multi-scale dictionaries. Fu et al. [32] utilized the multi-scale regions of an image to train a recurrent attention network for fine-grained recognition. Liu et al. [33] used multi-scale and multi-level network in a holistic manner to obtain hierarchical edge information. Similar to the inception network [34], Du et al. [35] fused different scale features from three varying filters.

The aforementioned CNN-based SR models build fine networks and have advanced the state-of-the-art performance on learning significant local detail information. The approach in [36] points out that too small receptive field resulting in the lack of enough global information to yield good visual results. To obtain fine local detail information, they often use small image patches for training (e.g., $33\times 33$ for SRCNN [17], $41\times 41$ for VDSR [19], $32\times 32$ for LapSRN [20], and $24\times 24$ for SRResnet/SRGAN [21]). A small receptive field only considers a limited range of information during SR tasks. This model lacks the capability of obtaining global and contextual information for SR. On the contrary, Zeiler et al. [37] visualized convolutional network to indicate that different network layers have varying roles in representing the features that simulate the ventral pathway to enhance their performance [38,39,40,41,42]. They indicated that hierarchical features of different scales effectively improve the capability of acquiring global information.

Inspired by the observation of “look closer to see better” [32], we propose a flexible and versatile multi-scale residual deep neural network for satellite image SR, named MRNN, for the hierarchical reconstruction of satellite imagery with HR detail information. In this network, multi-scale receptive fields are similar to the observation from different distances by human eyes. We extract three scales of the image at large- (large-kernel-size network, for global information), middle- (middle-kernel-size network, for contextual information) and small-scale (small-kernel-size network, for fine local information) features to represent the multi-scale information of images. In comparison with traditional neural networks, MRNN fuses the residual information rather than intermediate features. Thus, the fusion network fuses all scales of residual information to improve the high-frequency details.

The contributions of this study are highlighted as follows: (i) The use of MRNN is proposed for satellite image SR. The proposed network contains three parts, namely multi-scale feature extraction; parallel small-, middle-, and large-scale; and residual fusion networks. The proposed multi-scale neural network leverages SR performance on the basis of “look in multi-scale to see better”. (ii) The proposed residual enhancement and fusion networks effectively enhance the high-frequency information of satellite images in SR tasks. The fusion network refines fine edge/detail textures, thereby improving the details of the satellite image.

The remainder of this paper is organized as follows. In Section 2, we describe the framework of the proposed method. In Section 3, comparison is presented among the proposed method and some representative SR methods. The discussion and conclusion of this study are given in Section 4 and Section 5, respectively.

2. Satellite Imagery SR Based on Multi-Scale Residual Neural Network

We use image saliency to show the difference among various image sizes and emphasize the role of multi-scale images. Image saliency [43] is an important visual feature in an image and emphasizes the importance degree of a region for human eye perception. The brightness of a saliency map represents the importance of object parts. The saliency map S is formulated as:

$$S\left(x,y\right)={∥{\mathit{I}}_{\mu }-{\mathit{I}}_{\omega hc}\left(x,y\right)∥}_2^2,$$

(1)

where ${\mathit{I}}_{\mu }$ is the mean image feature vector, and ${\mathit{I}}_{\omega hc}\left(x,y\right)$ is the corresponding image pixel vector value at position $(x,y)$ in the Gaussian blurred version (using a $5\times 5$ separable binomial kernel) of the original image.

Figure 1 displays three sizes of image patches, namely large- ($91\times 91$), middle- ($61\times 61$), and small-sized ($41\times 41$) image patches. In the $91\times 91$ image patch, the saliency map focuses on the global information in the image, such as the outline of a building. For the $61\times 61$ image patch, which contains further contextual information, the saliency map focuses on building parts and street lines. For the $41\times 41$ image patch, which has a small receptive field, only local information is observed, and global information is neglected. Here, long-distance observation experience can be reviewed; global configuration information, such as position and outward appearance, can be observed when we are far from the observed objects, and no detailed information is included. For additional details, we focus on local information, such as the decoration and color of a building, as we approach. This observation is a good illustration of the role of multi-scale information in visual observation. Therefore, image reconstruction on only single-scale image patches cannot simultaneously and effectively recover the global and local information of the object.

We propose a novel multi-scale residual network, whose structure is shown in Figure 2. We establish three adaptive networks with different scale features to predict their high-frequency residual information in different scales for satellite images. Thus, we use residual images with varying scales to merge their high-frequency by utilizing a residual fusion network. As the pixel value in the residual image is small, we use the ImageEnhance module of Python Imaging Library https://github.com/python-pillow/Pillow to conduct enhanced contrast processing of images. The enhanced image $blend\_img$ is given by

$$blend\_img=img1\times \left(1-\lambda \right)+img2\times \lambda ,$$

(2)

where $img1$ is the original image, and the enhancement factor $\lambda =10$ represents the weight of the image blend. $img2$ is a generated image, whose pixel value is 0.5 plus the average value of $img1$. The greater is the $\lambda$, the greater is the contrast of the image.

For a pair of training datasets ${\{{\overline{\mathit{X}}}_i,{\mathit{Y}}_i\}}_{i=1}^M$, where LR image ${\overline{\mathit{X}}}_i\in {\Re }^{h\times w}$ and HR image ${\mathit{Y}}_i\in {\Re }^{ht\times wt}$, t is the amplification factor, i denotes the sample index, and M refers to the number of training samples. The LR image ${\overline{\mathit{X}}}_i\in {\Re }^{h\times w}$ is interpolated to the HR image size with bicubic kernel as ${\mathit{X}}_i\in {\Re }^{ht\times wt}$, and the tensor version of the training dataset is rewritten as ${\{{\mathit{x}}_i,{\mathit{y}}_i\}}_{i=1}^M$. The superscript represents the type of network, and the subscript indicates the number of layers. Superscripts K3, K5, K7, C, and F represent the K3-network, K5-network, K7-network, Concat operation, and residual fusion network, respectively. The sampling of patches with different sizes results in various numbers of patches in each scale. However, all training sample sets share the same training set ${\{{\mathit{x}}_i,{\mathit{y}}_i\}}_{i=1}^M$. The number of image patches is calculated as follows:

$$N_{S_D}=⌊ht/S_D⌋\ast ⌊wt/S_D⌋\ast M,$$

(3)

where $\lfloor \rfloor$ is the floor function and $S_D$ indicates the size of receptive field of the D-layer network. Image patches $41\times 41$ and $61\times 61$ are acquired on the basis of the center point of the $91\times 91$ image patch (for additional details, see Figure 2). LR and HR image patch pairs with different scales are defined as ${\left\{{\mathit{x}}_j^{K3},{\mathit{y}}_j^{K3}\right\}}_{j=1}^{N_{41}}$, ${\left\{{\mathit{x}}_j^{K5},{\mathit{y}}_j^{K5}\right\}}_{j=1}^{N_{61}}$, and ${\left\{{\mathit{x}}_j^{K7},{\mathit{y}}_j^{K7}\right\}}_{j=1}^{N_{91}}$, which have patch sizes of $41\times 41$, $61\times 61$ and $91\times 91$ pixels, respectively. j is the index of the image patches, and $N_{41}$, $N_{61}$, and $N_{91}$ denote the numbers of patches. Considering residual fusion, we use the patch center point to anchor three different size patches; thus, $N_{41}=N_{61}=N_{91}$.

2.1. Multi-Scale SR

We use three different scales of networks to simulate SR with different depths. The network depths are $D_{k3}$, $D_{k5}$, and $D_{k7}$. Parameter D is fine tuned according to the method in Section 3. In the K3-network, the convolution filter is defined as $k=3$. The residual map of the K3-network at the patch level is defined as follows:

$$f^{K3}\left({\mathit{x}}_j^{K3}\right)={\mathit{W}}_{20}^{K3}\ast H_{19}^{K3}\left({\mathit{x}}_j^{K3}\right)+b_{20}^{K3},$$

(4)

where $f^{K3}\left({\mathit{x}}_j^{K3}\right)$ is the predicted residual patch with size $41\times 41$; ${\mathit{W}}_{20}^{K3}$ indicates the weight matrix with size $64\times 3\times 3\times 1$; ${\mathit{b}}_{20}^{K3}$ denotes the bias with size $1\times 1$; $H_{19}^{K3}$ represents the generated feature maps of the 19th layers by an activation ReLU, which is composed of 64 feature maps; and j refers to the index of image patches.

For the K5-network, the size of its convolution kernel is $5\times 5$ pixels. For K7-network, the filter kernel size is $7\times 7$ pixels. We use the same method to calculate the size of the input image patch. Their residual maps are calculated as follows:

$$f^{K5}\left({\mathit{x}}_j^{K5}\right)={\mathit{W}}_{15}^{K5}\ast H_{14}^{K5}\left({\mathit{x}}_j^{K5}\right)+b_{15}^{K5},$$

(5)

$$f^{K7}\left({\mathit{x}}_j^{K7}\right)={\mathit{W}}_{15}^{K7}\ast H_{14}^{K7}\left({\mathit{x}}_j^{K7}\right)+b_{15}^{K7},$$

(6)

where ${\mathit{W}}_{15}^{K5}$ has a size of $64\times 5\times 5\times 1$; ${\mathit{W}}_{15}^{K7}$ has a size of $64\times 7\times 7\times 1$; the size of $b_{15}^{K7}$ and $b_{15}^{K5}$ is $1\times 1$; and j denotes the index of the image patches. $H_{14}^{K5}$ and $H_{14}^{K7}$ represent the feature maps of the 14th layers by the K5- and K7-networks, respectively.

2.2. Residual Fusion Network

To realize the complementarity of different scales of information, the global information of an object is described by large-scale information, and, the closer you look, the better the hierarchical details become. We use a fusion network for multi-scale residual fusion.

$$f^C\left(\mathit{x}\right)=Concat(f^{K3}{\left({\mathit{x}}_j^{K3}\right)}_r,f^{K5}{\left({\mathit{x}}_j^{K5}\right)}_r,f^{K7}{\left({\mathit{x}}_j^{K7}\right)}_r)=\left[f^{K3}\right({\mathit{x}}_j^{K3}{)}_r,f^{K5}\left({\mathit{x}}_j^{K5}{)}_r,f^{K7}\right({\mathit{x}}_j^{K7}{)}_r],$$

(7)

where $f^{K3}{\left({\mathit{x}}_j^{K3}\right)}_r$, $f^{K5}{\left({\mathit{x}}_j^{K5}\right)}_r$, and $f^{K7}{\left({\mathit{x}}_j^{K7}\right)}_r$ are the residual maps with the removal of border from the outputs of the three differently scaled networks. $f^C\left(\mathit{x}\right)$ represents the combined three layers of residual maps. The Concat function cascades the multi-scale residual maps in the third dimension (connect three tensors). Regardless of the same input $\mathit{x}$, the outputs of K3-, K5-, and K7-networks are different because they reconstruct their residual information through their own scales. To fuse different scales of residual information, we use a simple two-layer network to fuse three channel information. A $1\times 1$ convolution kernel is a linear combination of each pixel on different channels. The $1\times 1$ convolution kernel is used to fuse the residual feature maps. The cross-channel information interaction among different scales of information is consistent with the hierarchical visual cognition mechanism. We can obtain the final fusion residual as follows:

$$R^F\left(f^C\left(\mathit{x}\right)\right)={\mathit{W}}_2^F\ast H_1^F\left(f^C\left(\mathit{x}\right)\right)+b_2^F,$$

(8)

where ${\mathit{W}}_2^F$ is the second layer weight matrix, $b_2^F$ represents its bias, $f^C\left(\mathit{x}\right)$ denotes the input multi-scale residual maps, and $R^F\left(\mathit{x}\right)$ indicates the final fused output residual map. Thus, the final HR image $\widehat{y}$ is as follows:

$$\widehat{y}=R^F\left(f^C\left(\mathit{x}\right)\right)+\mathit{x}.$$

(9)

2.3. Loss Function

We define the loss function with mean squared error (MSE) as the objective function. In MRNN, we formulate the overall loss function as follows:

$$\begin{array}{c}\begin{array}{c}\mathrm{Loss}=\alpha \sum _{j=1}^{N_{91}}{∥{\mathit{y}}_j^{K3}-{\mathit{x}}_j^{K3}-f^{K3}\left({\mathit{x}}_j^{K3}\right)∥}_2^2+\beta \sum _{j=1}^{N_{91}}∥{\mathit{y}}_j^{K5}\hfill \\ {-{\mathit{x}}_j^{K5}-f^{K5}\left({\mathit{x}}_j^{K5}\right)∥}_2^2+\chi \sum _{j=1}^{N_{91}}{∥{\mathit{y}}_j^{K7}-{\mathit{x}}_j^{K7}-f^{K7}\left({\mathit{x}}_j^{K7}\right)∥}_2^2\hfill \end{array}\hfill \\ {+\delta ∥{\mathit{y}}_j^{K3}-{\mathit{x}}_j^{K3}-R^F\left(f^C\left(\mathit{x}\right)\right)∥}_2^2,\hfill \end{array}$$

(10)

where the first three terms are the losses of the multi-scale residual networks (K3-, K5-, and K7-networks). The last term represents the residual fusion loss. We simply set $\alpha =\beta =\chi =\delta =1$. We use a two-step method to train the network. Initially, we parallel-train three SR networks with differently-scaled patches. Then, we determine the fusion loss for the second time on the basis of the contacted residual maps.

A gradient descent method is used to optimize the network parameters by back propagation. Convolution operations reduce the size of the feature map. We maintain many edge pixels by padding zero to infer the center pixel accurately and ensure that all feature maps have the same size to preserve the information on the edge of the image patch.

3. Experiments

3.1. Experimental Data

The learning-based super-resolution methods learn the missing high-frequency information of LR images from the prior information provided in the training data. Generally, the more training data there are, the better reconstruction effect can be obtained by SR methods. In addition, the performance of the SR reconstruction method is also related to the similarity of the test image to the training image. If the test image is close to the statistical characteristics of the training images, it is more likely to get a good reconstruction result. At this point, there may be fewer training samples to get good results. On the contrary, when the statistical characteristics of the test image and the training image are greatly different, it is difficult to achieve a satisfactory result even using a large-scale training set. To verify the performance of MRNN, we conducted experiments on two satellite image datasets, namely, SpaceNet image and NWPU-RESISC45, to ensure that all algorithms used the same amount of training data. The SpaceNet satellite image dataset https://spacenetchallenge.github.io/AOI_Lists/AOI_1_Rio.html includes five areas in Rio de Janeiro, Paris, Las Vegas, Shanghai, and Khartoum, which are collected from DigitalGlobe’s WorldView-2 satellite and published publicly at Amazon. The complete satellite image of Rio de Janeiro (the spatial resolution is 0.5 m) has the highest resolution image with 2.8 M × 2.6 M pixels, and is divided into 6540 non-overlapping HR image patches with $436\times 404$ pixels, and the main contents of interest in the image are buildings and roads. In total, 2080 images of buildings were randomly selected from these image patches, of which 2000, 40, and 40 images were used as the training set, validation set, and test samples, respectively.

The NWPU-RESISC45 dataset http://pan.baidu.com/s/1mifR6tU [44] is a publicly available benchmark for remote sensing image scene classification (RESISC), created by Northwestern Polytechnical University (NWPU). This dataset covers 45 classes with 700 images in each class. We randomly selected 52 images from each class, of which 50 were used for training and the rest for testing. The HR image size is $256\times 256$ pixels. The spatial resolution of NWPU-RESISC45 varies from approximately 30 m to 0.2 m [44]. Images in the NWPU-RESISC45 dataset, compared with the SpaceNet dataset, have complicated and erratic imaging conditions, including various weather, seasons, and lighting conditions. These factors pose a huge difficulty for SR methods.

Image degradation is a very complex process to be modeled by some filter and down-sampling operators. Here, we interpolated the HR image with bicubic kernel into its LR version with scaling factor t. In the current works (for example, all the comparison methods in our work [17,20,21,23,45]), the most commonly used image degradation is the bicubic downsampling. Since learning-based super-resolution algorithms learn the mapping relationship between low-resolution and high-resolution images, the bicubic degradation is the fairest approach for comparison. Complex imaging degradation model will be investigated in future research. In the testing process, the images did not need to be partitioned.

Peak signal to noise ratio (PSNR) and structural similarity (SSIM) [46] (with default parameters) describe the similarity between the reconstructed and original images in terms of the image. Recent studies [47] have shown that feature similarity (FSIM) [47] and visual information fidelity (VIF) [48] are further consistent with the subjective results. Rectangular-normalized superpixel entropy index (RSEI) [49] (with default parameters) https://github.com/jiaming-wang/RSEI obtains further accurate image evaluation results by introducing the spatial structure of the image. Mutual information (MI) can express the dependence degree of the information between the images in terms of information. The higher is the MI score, the more substantial is the dependence and the higher is the similarity between images. The mutual information between patches $\mathit{y}$ and $\widehat{\mathit{y}}$ is defined as follows:

$$MI(\widehat{\mathit{y}};\mathit{y})=\sum _{q\in \mathit{y}}\sum _{g\in \widehat{\mathit{y}}}P(q,g)log\frac{P(q,g)}{P\left(q\right)P\left(g\right)},$$

(11)

where q and g represent the gray-scale values, $P\left(g\right)$ denotes the ratio of the number of pixels of the gray value that is g to the increased image, and $P(q,g)$ is the joint distribution function of q and g.

We define the information gain between SR image $\widehat{\mathit{y}}$ and LR image $\mathit{x}$ relative to HR image $\mathit{y}$ as follows:

$$GMI(\widehat{\mathit{y}};\mathit{y})=\frac{MI(\widehat{\mathit{y}};\mathit{y})}{MI(\mathit{x};\mathit{y})}.$$

(12)

All image quality assessment metrics only consider the Y component of the YCbCr color space.

3.2. Training Parameters

The proposed network is an end-to-end network, where each sub-network must train for 80 epochs as the pre-training network. The entire network is trained for 10 epochs.

Considering the deep network layer, the algorithm uses learning rate attenuation. We followed Kim et al. [19] for setting hyper-parameters: the learning rate was initialized to 0.1, the learning rate decreased by 1/10 every 20 epochs, and the network’s momentum was 0.9. To avoid over-fitting, we used regularized ${\ell }_2$-norm, and its weight decay was 0.0001. For the K3-residual learning network, we set the step size to 1 with a padding size of 1. For the K5-network, the step size was equal to 1 with a padding size of 2. For the K7-network, we set the step size to 1 and padding size to 3. We applied the MSRA method [50] to initialize the weights, that is, satisfying the Gaussian distribution whose mean value is 0, utilizing a variance of $\sqrt{\frac{2}{n}}$ (n is the batch size), and a constant to initialize the bias term with initial value 0. We initially converted the RGB image to the YCbCr color space and then reconstructed the Y channel. After the reconstruction, the Y channel image was restored to the RGB color space. We implemented the MRNN model using the Caffe library [51]. Training the MRNN roughly took 10 h with four 1080Ti GPUs.

3.3. Complementarity Analysis of Multi-Scale Residual

If there were less overlap between different scale residual information, it would mean that the complementarity of residual information between different scales is better [52]. Therefore, in this section, we show the distributions of residual information on different scales. We selected 15 representative LR images $\left\{{\mathit{x}}_i\right\}(i=1,\dots ,15)$ and corresponding HR image $\left\{{\mathit{y}}_i\right\}$ from SpaceNet image datasets with the same configuration of Section 3.6. The reconstruction residual maps of multi-scale networks are $\left\{f^j\left({\mathit{x}}_i\right)\right\}(j=K3,K5,K7)$ for a total of 45 residual images. The estimation residual error map was defined as ${\mathit{erm}}_i^j=f^j\left({\mathit{x}}_i\right)-({\mathit{y}}_i-{\mathit{x}}_i)(i=1,\dots ,15andj=K3,K5,K7)$, and we projected them into 3D and 2D residual feature spaces through principal component analysis (PCA), as shown in Figure 3. The distribution maps of 2D and 3D feature space show that multi-scale networks provide different estimation residual errors. This observation also proved that they are complementary. The overlap observed in Figure 3B covers a sufficiently large feature space, even if only three parallel networks are used. Therefore, additional parallel networks would only increase overlap.

The distribution maps in Figure 3 cannot clearly describe the complementary patterns of multi-scale residual. Therefore, we implemented the clustering of data by k-means and obtained their distribution of 2D feature space, as shown in Figure 4. We name the four patterns as “$s+m+l$”, “$s+m/l$”, “$m+s/l$”, and “$l+s/m$”. Pattern “$s+m+l$” represents the best case, that is, the high-frequency information of three scales is complementary between any two. The latter three patterns can be classified as: the information of two scales is considerably common, but a complementary relationship also exists, whereas the other scale complements them. This behavior effectively demonstrates the complementarity between multi-scale residuals.

We performed quantitative validation as follows:

$$\begin{array}{c}{C}_{{erm}^j}=card(\left|{\mathit{erm}}^j\right|>t),\hfill \\ C_{overlap}=card(\left|{\mathit{erm}}^{K3}\right|\&\left|{\mathit{erm}}^{K5}\right|\&\left|{\mathit{erm}}^{K7}\right|>t),\hfill \\ wherej=K3,K5,K7,\hfill \end{array}$$

(13)

where $abs(.)$ represents the absolute value of the matrix in an element-wise manner. The function $card(.)$ can count the number of nonzero elements in a matrix. $C_{{erm}^j}$ represents the number of elements whose values are greater than threshold t. $C_{overlap}$ denotes the number of above elements at the same locations in three error residual maps. We refer to Wang et al. [52] and set $t=9$ to represent high-value components (high-frequency information signals). Figure 5 plots the bar. The blue bar represents the error only from the K3-network, and the green and red bars indicate the errors only from the K5- and K7-networks, respectively. $C_{overlap}$ is the purple bar. $\forall j\in \left\{\mathrm{K}3,\mathrm{k}5,\mathrm{k}7\right\},C_{overlap}<C_{{erm}^j}$, and networks of different scales play different roles in the proposed method.

Statistical data and qualitative assessments prove that high-frequency information learned by multi-scale networks is complementary. This case is the reason we fuse multi-scale residual maps for improving reconstruction performance.

3.4. Performance and Model Trade-Offs

We configured the multi-scale residual network to different depths and compared their performance. We set D at 5, 10, 15, 20, and 25 to test the network performance. The input image patch size changed when the network depth changed. We used PSNR to measure the network performance, as shown in Figure 6. For the K3-network, the performance was optimal when D was 20. For K5- and K7- networks, the performance of networks was optimal when D was 15. The receptive field $S_D\times S_D$ of the D-layer network is defined as $S_D=(k-1)\times D+1$, and k is the kernel size.

3.5. Visualizing the Learned Filters and Feature Maps

The experiments presented in the previous section showed that three different depths of $3\times 3$ networks can replace MRNN. The results prove that “deeper is not better” in certain low-level vision tasks. We would like K5- and K7-networks to learn contextual and global information to compensate for the lack of information in the K3-network. Therefore, we visualize the networks to consider the role of differently-scaled networks in this section.

In the recognition task, the features learned by the network exhibit hierarchical features. Deep features are more discriminative than shallow features, such as color and edge. Therefore, horizontal visualization is suitable for describing the recognition process from low to high level. The image restoration is different from the recognition task. To explore the role of differently scaled networks, we longitudinally visualize the MRNN, that is, the filters and feature maps of the penultimate layer of the differently scaled networks.

A large difference is observed in the complexity of patterns from the filters. Figure 7 represents the feature maps. The larger the filters are, the less local detail information is represented in the feature maps. The smaller are the filters, the more apparent is the detail information in the feature maps.

Overall, we observe that differently scaled networks have their own advantages on various scale objects. For example, a large-scale network performs efficiently on global configuration, a middle-scale network is good at contextual information, and a small-scale network performs well in local detail information. K3-, K5-, and K7-networks have different levels of functionality in the network. A single-scale network cannot simultaneously learn different scales of information. Thus, the multi-scale information should be fused to improve image reconstruction performance.

3.6. Performance Comparison with State-Of-The-Art SR Algorithms

We conducted subjective qualitative and quantitative analyses on the reconstructed images by using PSNR, SSIM, FSIM, VIF, RSEI, and GMI. To verify the effectiveness of our algorithm, we compared MRNN with the following state-of-the-art SR algorithms:

• SelfExSR [45] is the best performing algorithm based on self-similarity based SR.
• SRCNN [17] is a classic deep-learning based approach, which first uses CNN for SR task.
• LapSRN [20] is the most famous multi-scale SR algorithm based on deep learning.
• VISR [23] is the best performing of satellite image SR algorithm via CNN.
• SRResnet [21] is an excellent depth network algorithm with high computing efficiency and high visual fidelity.

These algorithms were implemented using their public source codes and available parameters provided by the authors, and all images were down-sampled by using the same bicubic kernel of MATLAB. For a fair comparison, we trained all these algorithms with the same database configuration and evaluated the same satellite images with the proposed network.

Figure 8 shows the PSNR, SSIM, FSIM, VIF, RSEI, and GMI of all 40 testing images. MRNN obtained improved reconstruction results. The corresponding significance levels were 100%, 100%, 97.5%, 100%, 100%, and 95%, respectively. The difference in score between MRNN and other methods was statistically significant.

Table 1. Average results of PSNR, SSIM, MI, and GMI on the SpaceNet dataset with scale factor of 4. Bold indicates the best performance.

Eval. Mat	Bicubic	SelfExSR	SRCNN	LapSRN	VISR	SRResnet	MRNN
PSNR	25.13	26.20	26.53	26.52	26.74	26.38	27.02
SSIM	0.7262	0.7613	0.7675	0.7672	0.7793	0.7738	0.7894
FSIM	0.9097	0.9448	0.9501	0.9503	0.9570	0.9546	0.9575
VIF	0.3272	0.3848	0.3917	0.3864	0.4021	0.3964	0.4124
RSEI	0.3590	0.3760	0.3782	0.3785	0.3828	0.3797	0.3859
MI	5.0582	5.1229	5.0980	5.1040	5.1114	5.0972	5.1424
GMI	1.0000	1.0128	1.0079	1.0091	1.0106	1.0078	1.0167

and

Table 2. Average results of PSNR, SSIM, MI, and GMI on the NWPU-RESISC45 dataset with scale factor of 4. Bold indicates the best performance.

Eval. Mat	Bicubic	SelfExSR	SRCNN	LapSRN	VISR	SRResnet	MRNN
PSNR	25.57	28.48	28.49	28.81	28.77	28.80	28.93
SSIM	0.6920	0.7403	0.7378	0.7578	0.7503	0.7564	0.7580
FSIM	0.7872	0.8384	0.8336	0.8400	0.8386	0.8419	0.8461
VIF	0.2968	0.3499	0.3464	0.3612	0.3661	0.3678	0.3685
RSEI	0.3572	0.3740	0.3760	0.3766	0.3777	0.3807	0.3837
MI	4.4972	4.5807	4.5507	4.5517	4.5478	4.5808	4.6210
GMI	1.0000	1.0191	1.0122	1.0125	1.0117	1.0191	1.0283

show a considerable quantitative advantage of the proposed method compared with cutting-edge deep learning based algorithms. This finding indicates that residual multi-scale networks are relatively effective in learning different scales of content and structure, and they restore image information effectiveness by using a deeper and flatter network than those used by competing algorithms.

For simple observation, we amplified the representative scale object in randomly selected reconstructed image for comparison. As shown in Figure 9, we selected a roof (small-scale object), building (middle-scale object), and street corner (large-scale object) to show the SR performance. For the examples shown in Figure 9, our method produced sharper edges and finer details than the other methods for all object scales. In addition, our method produced sharper edges and finer details than LapSRN for all object scales. This condition confirms that MRNN fuses multi-scale residual information to enhance visual performance. Figure 10 shows a further intuitive result that only our method can restore a clear outline.

3.7. Time Complexity

Figure 11 shows the running time of all algorithms. The running time of the traditional algorithm is longer than that of deep learning algorithms and has no training phase. MRNN is a parallel network with three different scales and does not increase the time complexity of the network, especially when the network is complex. Although LapSRN has a better running time performance, its PSNR is lower than that of MRNN. Our method is slightly slower than VISR in terms of running time. However, MRNN has improved PSNR, SSIM, FSIM, VIF, RSEI, and GMI. We implemented all algorithms in the experiments under the same hardware configuration: Intel Core i7-6700 K CPU @4.00 GHz, NVIDIA GTX1080 8 GB RAM.

4. Discussions

4.1. Multi-Scale Prior Information

Lai et al. [20] proposed a progressive SR method to super-resolve images gradually. A Laplacian pyramid is used in the generative network for SR. A residual recurrent network is adopted to predict the output information in each pyramid level. Here, LapSRN designs a multi-scale training strategy, which trains multi-scale combinations as 2×, 4×, and 8× in one net. This process involves the addition of multi-scale training pairs to cover different scale samples. Many differences are observed between LapSRN and MRNN. First, LapSRN directly performs multi-scale information fusion in the feature domain, whereas MRNN constructs multi-scale parallel networks and performs multi-scale information fusion in the residual domain. The residual map is closely related to the high frequency information of the image, which is the purpose of SR. Second, LapSRN can perform SR tasks at different scales of factor one-shot, but it ignores the multi-scale information in the input image. MRNN completely investigates the multi-scale information of the input image in SR at fixed-scale factors. On the basis of the experimental results, the fusion of multi-scale residual information has a better performance than LapSRN at a scale factor of 4.

4.2. Residual Learning Versus Pixel Learning

VISR [23] uses a self-similar padding instead of zero padding to avoid the addition of unnecessary information. Therefore, VISR is performed in the pixel domain, such as SRCNN. The reconstruction in the pixel domain focuses on the low- and middle-frequency information in the image. However, SR infers the missing high-frequency information. Furthermore, the residual values of images are frequently small or zero, and the residual network has consistently less calculation burden than pixel learning. The recovery of high-frequency information on satellite images can improve recognition performance. At this point, residual learning is further suitable for satellite image SR scenarios. In addition, the experimental results confirm this inference in terms of subjective and objective image qualities.

4.3. Subpixel Network Versus Pixel Network

SRResnet [21] directly divides LR images into small image patches. Similar to LapSRN and ESPCNN [18], these networks directly use LR inputs (subpixels) to learn mapping functions. Subpixel networks simulate the degradation process and are more efficient than pixel-based networks. By contrast, pixel networks interpolate LR inputs into the same size of HR samples and use the residual information in networks. Residual recursive networks are assumed to overcome the vanishing problem for improving network performance. Thus, subpixel and pixel networks have their own advantages. On the basis of the experimental data in the SpaceNet database, the pixel network outperforms its subpixel competitors.

4.4. Applicability of the Proposed Method

We conducted experiments on the Jilin-1 satellite image to further illustrate the applicability of the proposed algorithm. The imaging environment and resolution of the test image are different from those of the training datasets (NWPU-RESISC45). The size of LR Jilin-1 satellite image is $408\times 204$ pixels. Figure 12 shows the reconstruction results obtained from our proposed approaches and the comparison methods. Considering the absence of ground truth image, we introduce mean gradient (MG) to calculate the sharpness of the SR image. MG is defined as follows:

$$\mathrm{MG}=\left|gr{d}_x\left(\mathit{y}\right)\right|+\left|gr{d}_y\left(\mathit{y}\right)\right|,$$

(14)

where $gr{d}_x\left(\mathit{y}\right)$ and $gr{d}_y\left(\mathit{y}\right)$ are the gradients of image $\mathit{y}$ on the x- and the y-axes, respectively. The proposed MRNN recovers sharp edges, and enjoys the first MG scores. The comparison results of real video satellite images show the applicability of the proposed method. Considering the image characteristics between satellites, we introduce GAN to learn the cross-domain degradation model for solving the real-world SR problems in the future.

5. Conclusions

This paper presents a multi-scale residual CNN, namely MRNN, based on the characteristics of satellite images, for enhancing SR performance. It first extracts different sizes of patches from LR satellite images. Then, multi-scale deep residual neural networks are applied to simulate differently sized receptive fields for acquiring different levels of information. Then, a fusion network is used to refine the multi-scale features. Based on the proposed novel network, reasonably accurate high-frequency information, such as edges and textures, can be obtained by complementing the residual information at different scales. The experimental results on the SpaceNet database show that the proposed MRNN effectively enhanced the high-frequency information in the reconstructed images. MRNN also exhibited better subjective and objective image qualities than several state-of-the-art deep-learning-based SR algorithms for satellite images. MRNN is mainly designed for true color satellite images SR. As is known, multi-spectral images have higher spectral resolution but lower spatial resolution. It would be very interesting to investigate the fusion of the multi-spectral images and the true color images in the proposed framework to improve the visualized quality of the multi-spectral images in the future.