Multi-Frame Super-Resolution of Gaofen-4 Remote Sensing Images

Gaofen-4 is China’s first geosynchronous orbit high-definition optical imaging satellite with extremely high temporal resolution. The features of staring imaging and high temporal resolution enable the super-resolution of multiple images of the same scene. In this paper, we propose a super-resolution (SR) technique to reconstruct a higher-resolution image from multiple low-resolution (LR) satellite images. The method first performs image registration in both the spatial and range domains. Then the point spread function (PSF) of LR images is parameterized by a Gaussian function and estimated by a blind deconvolution algorithm based on the maximum a posteriori (MAP). Finally, the high-resolution (HR) image is reconstructed by a MAP-based SR algorithm. The MAP cost function includes a data fidelity term and a regularized term. The data fidelity term is in the L2 norm, and the regularized term employs the Huber-Markov prior which can reduce the noise and artifacts while preserving the image edges. Experiments with real Gaofen-4 images show that the reconstructed images are sharper and contain more details than Google Earth ones.


Introduction
Images with higher resolution are required in most electronic imaging applications, such as remote sensing, medical diagnostics, and video surveillance. During the past decades, considerable advances have been realized in imaging systems. However, the quality of images is still limited by the cost and the manufacturing technology [1]. Super-resolution (SR) is a promising digital image processing technique to construct high-resolution (HR) images from several blurred low-resolution (LR) images.
The basic idea of SR is that the LR images of the same scene contain different information because of relative sub-pixel shifts, thus, an HR image with higher spatial information can be reconstructed by image fusion. Sub-pixel motion can occur due to movement of local objects or vibrating of the imaging system, or even controlled micro-scanning [2,3]. Numerous SR algorithms have been proposed since the concept of SR was introduced by Tsai and Huang [4] in the frequency domain. However, most of researchers nowadays address the problem mainly in the spatial domain, because it is more flexible to model all kinds of image degradations [5]. Ur and Gross [6] proposed a non-uniform interpolation of multiple spatially-shifted LR images based on the generalized multichannel sampling theorem, followed by a deblurring process. Stark et al. [7] were the first to suggest an SR algorithm based on the projection onto convex sets approach (POCS). Patti et al. [8] extended the method to handle arbitrary sampling lattices and non-zero aperture time. The advantage of the POCS method is that it allows a convenient inclusion of prior knowledge into the reconstruction process. Irani and Peleg [9] estimated the HR image by back projecting the error between re-projected LR images and the observed LR images. For the ill-posed nature of the SR problem, regularization methods are widely used in SR reconstruction. Schultz and Stevenson [10] reconstructed a single HR video frame from a short image sequence using the maximum a posteriori (MAP) approach with a discontinuity-preserving image prior. They used a block-matching method to estimate the sub-pixel displacement vectors between LR images. A MAP framework for jointly estimating image registration parameters and the HR image was presented by Hardie et al. [11]. The sub-pixel displacement vectors were iteratively updated along with the HR image in a cyclic coordinate-descent optimization procedure. Zomet et al. [12] proposed a robust SR algorithm in which a median estimator was used as the gradient of the L 2 norm cost function. Farsiu et al. [13] used a cost function in the L 1 norm rather than modified the gradient directly to further increase the robustness. For remote sensing applications, researchers mainly focused on single-frame SR based on learning or sparse representation for lack of multiple images of the same scene [14][15][16][17]. Nevertheless, some researchers have made their contributions. Merino and Núñez [18] proposed a method named Super-Resolution Variable-Pixel Linear Reconstruction (SRVPLR) to reconstruct an HR image from multiple LR images acquired over a long period. This algorithm was derived from drizzle [19] which was designed to work with astronomical dithered under-sampled images. A prior histogram matching on LR images and high geometrical corrections were introduced to make it compatible with remotely-sensed images. Shen et al. [20] proposed an SR method to MODIS remote sensing images, where image registration in the range and spatial domains were performed before image reconstruction. Li et al. [21] proposed a MAP-based SR algorithm with a universal hidden Markov tree model, and tested it with Landsat7 panchromatic images captured on different days. Fan et al. [22] proposed a POCS-based SR algorithm with a point spread function (PSF) estimated by a slant knife-edge method. The method was tested on Airborne Digital Sensor 40 (ADS40) three-line array images and the overlapped areas of two adjacent Gaofen-2 images.
This paper presents a technique that performs multi-frame super-resolution of Gaofen-4 satellite images. Gaofen-4 is a staring imaging satellite with high temporal resolution, and its high orbit limits spatial resolution, therefore, multi-frame SR is possible and essential. First, we build the observation model for the Gaofen-4 satellite, then we propose the pipeline of our MAP-based SR algorithm. Our algorithm consists of three steps: image registration, PSF estimation, and image fusion. Image registration is performed both in the range and spatial domains, and the Gaussian PSFs of LR images are estimated via a modified MAP-based blind deconvolution algorithm. Finally, The HR image is estimated by an image fusion method based on the MAP method with a Huber image prior, while the registration parameters are further refined along with the HR image.
The rest of this paper is organized as follows: Section 2 addresses the SR problem and details our algorithm; Section 3 presents the experimental results with synthetic and real data; Section 4 concludes the paper.

Observation Model
To analyze the SR problem comprehensively, an image observation model must be formulated first. For notational simplicity, the image in our method is represented as a lexicographically-ordered vector. Let us consider the desired HR image o of size LM × LN, where L denotes the downsampling factor in the observation model, and M and N are the number of rows and columns of the LR images, respectively. Generally, the observation model-related HR image and the LR images is given by [1]: where M t is the subpixel motion matrix of size L 2 MN × L 2 MN; B t represents the blur matrix of size L 2 MN × L 2 MN; D t denotes the downsampling matrix of size MN × L 2 MN; n t is the lexicographically-ordered additive Gaussian noise vector with zero mean; and y t represents the lexicographically-ordered LR image vector of size MN × 1. The observation model represented by Equation (1) assumes an illumination invariant environment. However, this assumption may not be valid due to the variation of light conditions over time. The modified observation model is given by: where λ1, t and λ0, t denote the contrast gain and illumination offset, respectively, and I is a unit vector of size MN × 1. The motion matrix Mt could vary in time, while the downsampling matrix Dt and blur matrix Bt remain constant over time for most situations. Therefore, D and B are used instead of Dt and Bt, respectively. In the following of this section, we build the specific model for the Gaofen-4 satellite. Remote sensing sensors can be modeled by rational polynomial coefficients (RPC) geometrically. However, for registering images taken by the Gaofen-4 staring imaging satellite, the scene can be effectively modeled as a planar region for its high orbit when the region of interest (ROI) is not too large and, hence, a simple affine camera model is sufficient to describe the image geometry. The affine transformation contains six parameters and can be written as: where a 0 and b 0 denote translation along x and y directions respectively, and a 1 , a 2 , b 1 and b 2 are the parameters representing rotation, scale and deformation. Using matrix notation, the entire affine motion parameters are represented by θ = [a 0 , a 1 , We assume a shift-invariant, convolutional blur, which can be described by the PSF in the spatial domain or the modulation transfer function (MTF) in the frequency domain. The remote sensing images can be blurred by the atmosphere, the optical system, relative motion, and pixel integration. The overall PSF is the successive convolution of all of the individual terms: which respectively summarize the atmospheric, optical system, motion and sensor contributions. The atmosphere contributes through turbulence, light scattering, and attenuation. The spatial resolution of Gaofen-4 in the visible and near-infrared band is 50 m, and the contribution of turbulence at this spatial resolution is insignificant. The scattering and absorption of energy by the aerosols affects all spatial frequencies, therefore, causing edges in the image to be blurred. The aerosol MTF can be approximated by a Gaussian form exp −k(r/r c ) 2 for r below the cut-off frequency r c [23].
A limited aperture size of the optical system causes a diffraction blur, and the corresponding PSF is given by the Airy disk. Other optical aberrations, such as defocus and spherical aberrations, also contribute to the optical system, which can be seen as perturbations to the Airy disk. A general optical system can be expressed as a Gaussian function: where a and b denote the deviations of the two-dimensional Gaussian function along the x and y direction, respectively. The motion-related blur results from shift and vibrations. For pushbroom and whiskbroom sensors, the PSF can be characterized by a rectangular function along the scanning direction. We assume the motion blur of the Gaofen-4 satellite can be ignored for it is a staring imager equipped with focal-plane arrays. The finite size of the sensor pixels results in sensor blur which presents the integration of irradiance over pixels. The PSF is given by: where w is the size of a pixel. All PSF terms can are gathered into a single Gaussian term, which is given by: where we further assume the deviations along the x and y direction are equal and denoted by σ.

Image Registration
The success of the SR technique relies on the subpixel motion between images, and the image registration accuracy affects the quality of the reconstructed HR image significantly. Geometric registration and photometric registration are interactive. Most photometric registration approaches require accurate geometric registration, and geometric registration may fail when images are not aligned photometrically. We first perform geometric registration coarsely using speeded up robust features (SURF) [24], which is not invariant to image contrast gain and illumination bias. Then photometric and geometric registration are jointly performed using alternative minimization.
The relationship between the target image f T (x, y) and the reference image f R (x, y) is given by the following model: where x denotes the spatial coordinates; θ and λ are geometric and photometric parameters in matrix form, respectively; G and P are geometric and photometric models, respectively; and n(x) is the additive Gaussian noise. In this paper, we focus on the linear photometric model and affine motion. Thus, Equation (8) becomes: where λ 1 and λ 0 denote the contrast gain and illumination bias, respectively. The goal of the image registration is to estimate these parameters via minimization of an objective function in the least square sense: where g(x, y) is the transformation of the target image which is expressed as: and Ω represents the region of inliers. To detect inliers, we use a geometric criterion and a photometric criterion. Geometric criterion requires that the pixels of the warped image should have the same geometric field with the reference image, i.e., where M and N are the number of rows and columns of the reference image, respectively. The photometric criterion requires that the intensity difference between the reference image and the registered target image should be small enough, that is to say: where d is the intensity threshold and is set to 8. The solution of Equation (10) can be obtained by an optimization algorithm, such as the conjugate gradient descent algorithm or the Gaussian-Newton method. However, it is easy to fall into a local optimal solution when the misalignment between two images is large. Thus, SURF is employed to estimate the geometric parameters coarsely, then an iterative updating strategy is used to estimate the photometric and geometric parameters more precisely. Specifically, at each iteration, the affine motion parameters are firstly updated by a new variation of the Marquardt-Levenberg algorithm, which was proposed by Thévenaz et al. [25], given the photometric parameters are fixed, then the affine motion parameters are kept fixed, and the photometric parameters are updated by the conjugate gradient method. The gradients of F(θ, λ) with respect to λ 1 and λ 0 are given by: where f W (x, y) is the target image after geometric registration, but without photometric registration. The image registration algorithm flow is shown in Algorithm 1. In the feature-matching step, the random sample consensus (RANSAC) method [26] is used to eliminate matching outliers. The stop criterion of this image registration algorithm is that values of the estimated parameters do not change anymore, which can be expressed by: where ∑ means the sum of all elements;θ k represents the estimated geometric parameters after k times cycles, and so on. ε is the precision threshold; here, 10 −4 is used.
Algorithm 1. Image registration algorithm. Select one image y k as the reference image, and detect its SURF features. 4 for l = 1, 2, . . . , T(l = k), do 5 Detect SURF features of the LR image y l . 6 Match features of y l and y k , then estimate geometric parameters θ based on RANSAC. 7 do 8 Estimate photometric parameters λ using conjugate gradient method for fixed geometric parameters θ. 9 Estimate geometric parameters θ using the Marquardt-Levenberg method for fixed photometric parameters λ. 10 until λ and θ do not change anymore. 11 endfor

PSF Estimation
The PSF of remote sensing images can be measured after being launched, or estimated via sharp edge prediction. Sometimes, appropriate sharp edges do not exist, so the PSF estimation method should be employed. We propose a modified MAP-based multi-frame blind deconvolution algorithm to estimate the PSF, which was originally proposed by Matson et al. [27]. The forward model used in this algorithm is: additive Gaussian noise; and * denotes convolution. The cost function to generate estimates of the PSFs is given by: where γ is the regularized parameter and σ 2 n,t is the image noise variance of the t-th frame. The minimization of the cost function is performed via the conjugate gradient method in a cyclic coordinate-descent optimization procedure. Readers can refer to [27] for more details of this blind deconvolution algorithm.
To make the above algorithm compatible with our application, we use the bicubic upsampled registered LR images as the measurement images. In order to save computational resources, no more than four measurement images are used. The estimated PSFs are reparameterized by Gaussian PSFs, which are given by Equation (7). The mean deviation of these PSFs is used in image fusion. We use the mean deviation because we assume there are no significant differences between the PSFs.

Image Fusion
Mathematically, the SR problem is an inverse problem whose objective is to estimate an HR image from multiple observed blurry LR images. The MAP estimator is a powerful tool to solve this type of problem. The MAP optimization problem can be represented by: Under the assumption of additive Gaussian noise, the likelihood term p(y t |o) can be expressed as: The term p(o) represents the probability density of the HR image, which restricts the solution space and introduces prior knowledge to the reconstruction. In particular, we employ a Huber Markov random field (HMRF) as the image prior. HMRF encourages piecewise smoothness, and can preserve image edges as well. The expression of p(o) is given by: where S denotes a gradient operator; v denotes the prior strength; α is a parameter of the Huber function; and Z is a normalization factor. The Huber function ρ is defined as: where z i is the i-th element of the vector z. By substituting Equations (21)- (24) to Equation (20) and neglecting the constant term, the MAP estimator can be converted to a minimizing cost function, which is given by: The accuracy of the geometric and photometric registration is limited by the low quality of the LR images. Thus, an iterative update of the registration parameters, along with the HR image, is required. For the discrepancy of convergence rate of x, {θ t }, and {λ t }, a cyclic coordinate-descent optimization procedure is adopted. At each iteration, the HR image is updated by minimizing Equation (25) with respect to o using a scaled conjugate gradient optimization method, given the registration parameters are fixed; then the registration parameters are updated by minimizing the cost function with respect to the registration parameters, given the HR image estimate is fixed. This minimization is completed by an image registration process, which is the same with the registration method proposed in Section 2.2. However, feature-based registration is not needed anymore, and the current HR image estimate is used as the reference image instead. The algorithm stops when the value of the function does not change. The proposed algorithm as a whole system is illustrated in Figure 1. The algorithm consists of three stages: image registration, PSF estimation and image fusion.
The accuracy of the geometric and photometric registration is limited by the low quality of the LR images. Thus, an iterative update of the registration parameters, along with the HR image, is required. For the discrepancy of convergence rate of x , , a cyclic coordinate-descent optimization procedure is adopted. At each iteration, the HR image is updated by minimizing Equation (25) with respect to o using a scaled conjugate gradient optimization method, given the registration parameters are fixed; then the registration parameters are updated by minimizing the cost function with respect to the registration parameters, given the HR image estimate is fixed. This minimization is completed by an image registration process, which is the same with the registration method proposed in Section 2.2. However, feature-based registration is not needed anymore, and the current HR image estimate is used as the reference image instead. The algorithm stops when the value of the function does not change. The proposed algorithm as a whole system is illustrated in Figure 1. The algorithm consists of three stages: image registration, PSF estimation and image fusion.

Simulation Images
We first test our algorithm using synthetic images. The quantitative evaluation of the algorithm is achieved by computing the root mean square error (RMSE) and Structural Similarity (SSIM). RMSE is defined as:

Simulation Images
We first test our algorithm using synthetic images. The quantitative evaluation of the algorithm is achieved by computing the root mean square error (RMSE) and Structural Similarity (SSIM). RMSE is defined as: where o andô are the original image and the reconstructed HR image respectively, and m is the total number of pixels of the original image. The SSIM index is defined as: [28] SSI (v) Apply photometric parameters with uniformly-distributing contrast gain in the range of (0.9, 1.1) and bias in the range of (−0.1, 0.1) × max(o). (vi) Make the SNR of the images become 30 dB by adding white Gaussian noise with zero mean.
A QuickBird 2 remote sensing image with a size of 256 × 256 is used as the origin HR image, and the downsampling factor is equal to 2. We use six LR images to reconstruct an HR image. The Huber function parameter α is set to 0.04, and the prior strength parameter v is set to 10. The estimated deviation of the total PSF is 1.12. Our method is compared with the bicubic interpolation, the iterative back projection (IBP) approach proposed by Irani and Peleg [9] and Farsiu's robust SR algorithm [13]. Note that the original IBP and Farsiu's algorithm do not consider the variations of the illumination in the image registration stage, thus we use our image registration algorithm instead. Figure 2 shows the original image and the reconstructed images by different methods. Compared with the bicubic upsampled image, the HR images reconstructed by all SR methods are sharper and contain more details. The reconstructed images are further compared quantitatively in RMSE and SSIM, as is listed in Table 1. The performance of our method and IBP is close, and better than the bicubic interpolation and Farsiu's algorithm. Additionally, without photometric registration, many artifacts appear in the reconstructed HR image, as is shown in Figure 2d.

Real Gaofen-4 Images
The second experiment is carried out with ten 5-band Gaofen-4 images without ortho-rectification. The images were all captured on 3 March 2017, from 11:10:20 to 11:20:21. The original size of full-view images is 10,240 × 10,240, and we chip the region of interest (ROI) with a size of 256 × 256. This region is an offshore aquaculture area, and the images captured at 11:10:20 and 11:20:21 are shown in Figure 3.

Real Gaofen-4 Images
The second experiment is carried out with ten 5-band Gaofen-4 images without ortho-rectification. The images were all captured on 3 March 2017, from 11:10:20 to 11:20:21. The original size of full-view images is 10,240 × 10,240, and we chip the region of interest (ROI) with a size of 256 × 256. This region is an offshore aquaculture area, and the images captured at 11:10:20 and 11:20:21 are shown in Figure 3.  The image captured at 11:10:20 is used as the reference image, and the other images are aligned to it. Although the images were captured over a short period, photometric registration is needed. The criterion used to evaluate the registration accuracy is signal to noise ratio (SNR), which is defined as:

Real Gaofen-4 Images
The second experiment is carried out with ten 5-band Gaofen-4 images without ortho-rectification. The images were all captured on 3 March 2017, from 11:10:20 to 11:20:21. The original size of full-view images is 10,240 × 10,240, and we chip the region of interest (ROI) with a size of 256 × 256. This region is an offshore aquaculture area, and the images captured at 11:10:20 and 11:20:21 are shown in Figure 3. The image captured at 11:10:20 is used as the reference image, and the other images are aligned to it. Although the images were captured over a short period, photometric registration is needed. The criterion used to evaluate the registration accuracy is signal to noise ratio (SNR), which is defined as:  The image captured at 11:10:20 is used as the reference image, and the other images are aligned to it. Although the images were captured over a short period, photometric registration is needed. The criterion used to evaluate the registration accuracy is signal to noise ratio (SNR), which is defined as: where g(x) is the target image after registration, and f R (x) is the reference image. The registration results are listed in Table 2. The images are numbered according to their capture time and the first LR image is used as the reference image. All images have similar registration accuracy with both geometric and photometric registration while, without photometric registration, the SNRs are lower, and become lower and lower as time progresses. The Huber function parameter α is set to 0.04, and the prior strength parameter v is set to 10. The estimated deviation of the total PSF is 1.21. The experimental results are shown in Figure 4. with both geometric and photometric registration while, without photometric registration, the SNRs are lower, and become lower and lower as time progresses. The Huber function parameter α is set to 0.04, and the prior strength parameter v is set to 10. The estimated deviation of the total PSF is 1.21. The experimental results are shown in Figure 4.  To validate our reconstructed HR image, we chip an image of the same zone from Google Earth, which was captured on 16 September 2015. The reconstructed HR image is clearer and contains more details which corresponds to Google Earth. To see more clearly, the regions denoted by the rectangles in Figure 4 are shown in Figure 5 without zooming. It is clear that the HR image reconstructed by our method is the sharpest. To validate our reconstructed HR image, we chip an image of the same zone from Google Earth, which was captured on 16 September 2015. The reconstructed HR image is clearer and contains more details which corresponds to Google Earth. To see more clearly, the regions denoted by the rectangles in Figure 4 are shown in Figure 5 without zooming. It is clear that the HR image reconstructed by our method is the sharpest. Quantitatively evaluating the performance of the proposed algorithm is not easy in real image applications, for there are no true original images. A commonly used strategy to solve the problem is to treat the image of an inferior resolution level such as the panchromatic image as the true original image. But the strategy is not applicable for Gaofen-4 images, because the panchromatic and multi-spectral images of Gaofen-4 have the same spatial resolution. Theoretically, the reconstructed HR image should be sharper than the bicubic upsampled image. Thus we assume the upsampled image is the blurred version of the reconstructed HR image and estimate the PSF between them. The PSF estimation is completed by the PSF estimation algorithm proposed by Matson et al. which has been described in Section 2.3. Then the PSF is fitted to the Gaussian function. A larger deviation of the Gaussian function means a sharper reconstructed HR image. The estimated PSFs for IBP, Farsiu's algorithm and our MAP-based SR algorithm are shown in Figure 6. The deviations of the fitted PSFs are 1.0, 0.9 and 1.0 respectively, which means our method has a close performance with IBP, and is better than Farsiu's algorithm. We extract another hilly region of interest with a size of 256 × 256 to further test our algorithm. The parameters are set the same as the previous experiment. The estimated deviation of the total PSF is 1.19, which is almost same with the previous experiment. The experimental results are shown in Figure 7. The Google Earth image was captured on 20 April 2017. To see more clearly, the regions denoted by the rectangles in Figure 7 are shown in Figure 8. Compared with other methods, our method results in the sharpest edges. Quantitatively evaluating the performance of the proposed algorithm is not easy in real image applications, for there are no true original images. A commonly used strategy to solve the problem is to treat the image of an inferior resolution level such as the panchromatic image as the true original image. But the strategy is not applicable for Gaofen-4 images, because the panchromatic and multi-spectral images of Gaofen-4 have the same spatial resolution. Theoretically, the reconstructed HR image should be sharper than the bicubic upsampled image. Thus we assume the upsampled image is the blurred version of the reconstructed HR image and estimate the PSF between them. The PSF estimation is completed by the PSF estimation algorithm proposed by Matson et al. which has been described in Section 2.3. Then the PSF is fitted to the Gaussian function. A larger deviation of the Gaussian function means a sharper reconstructed HR image. The estimated PSFs for IBP, Farsiu's algorithm and our MAP-based SR algorithm are shown in Figure 6. The deviations of the fitted PSFs are 1.0, 0.9 and 1.0 respectively, which means our method has a close performance with IBP, and is better than Farsiu's algorithm. To validate our reconstructed HR image, we chip an image of the same zone from Google Earth, which was captured on 16 September 2015. The reconstructed HR image is clearer and contains more details which corresponds to Google Earth. To see more clearly, the regions denoted by the rectangles in Figure 4 are shown in Figure 5 without zooming. It is clear that the HR image reconstructed by our method is the sharpest. Quantitatively evaluating the performance of the proposed algorithm is not easy in real image applications, for there are no true original images. A commonly used strategy to solve the problem is to treat the image of an inferior resolution level such as the panchromatic image as the true original image. But the strategy is not applicable for Gaofen-4 images, because the panchromatic and multi-spectral images of Gaofen-4 have the same spatial resolution. Theoretically, the reconstructed HR image should be sharper than the bicubic upsampled image. Thus we assume the upsampled image is the blurred version of the reconstructed HR image and estimate the PSF between them. The PSF estimation is completed by the PSF estimation algorithm proposed by Matson et al. which has been described in Section 2.3. Then the PSF is fitted to the Gaussian function. A larger deviation of the Gaussian function means a sharper reconstructed HR image. The estimated PSFs for IBP, Farsiu's algorithm and our MAP-based SR algorithm are shown in Figure 6. The deviations of the fitted PSFs are 1.0, 0.9 and 1.0 respectively, which means our method has a close performance with IBP, and is better than Farsiu's algorithm. We extract another hilly region of interest with a size of 256 × 256 to further test our algorithm. The parameters are set the same as the previous experiment. The estimated deviation of the total PSF is 1.19, which is almost same with the previous experiment. The experimental results are shown in Figure 7. The Google Earth image was captured on 20 April 2017. To see more clearly, the regions denoted by the rectangles in Figure 7 are shown in Figure 8. Compared with other methods, our method results in the sharpest edges. We extract another hilly region of interest with a size of 256 × 256 to further test our algorithm. The parameters are set the same as the previous experiment. The estimated deviation of the total PSF is 1.19, which is almost same with the previous experiment. The experimental results are shown in Figure 7. The Google Earth image was captured on 20 April 2017. To see more clearly, the regions denoted by the rectangles in Figure 7 are shown in Figure 8. Compared with other methods, our method results in the sharpest edges.   The PSFs between the HR images and the upsampled reference image of the hilly area are shown in Figure 9. The deviations of the fitted PSFs for IBP, Farsiu's algorithm and our MAP-based SR algorithm are 0.98, 0.82 and 1.0 respectively, which means our method yields the best SR result.
The PSFs between the HR images and the upsampled reference image of the hilly area are shown in Figure 9. The deviations of the fitted PSFs for IBP, Farsiu's algorithm and our MAP-based SR algorithm are 0.98, 0.82 and 1.0 respectively, which means our method yields the best SR result.

Further Discussion
The deviation of the PSF between the HR image and the upsampled image is less than the deviation estimated in the PSF estimation stage, i.e., the blur is not removed completely. A visually better image can be generated by contrast manipulation plus sharpening after the SR reconstruction. We take the real Gaofen-4 images of the offshore aquaculture area as an example. The result of enhancing the HR image reconstructed by our SR method (Figure 4c) is shown in Figure 10a. For comparison, the reference LR image is enhanced directly without SR reconstruction, and the result is shown in Figure 10b. The regions denoted by the rectangles in Figure 10a,b are shown in Figure  10c,d. Comparing Figures 4, 5 and 10, we note that the enhanced LR image is visually better than the LR image and the reconstructed HR image, but contains less details. That is to say, the SR algorithm reveals more information contained in the LR images.

Conclusions
In this paper, we propose an SR algorithm to reconstruct an HR image from multiple LR images captured by the Gaofen-4 staring imaging satellite. An appropriate observation model with photometric parameters is built. The registration with photometric adjustment improves the registration accuracy and the quality of the final reconstructed HR image. A MAP-based blind deconvolution algorithm is modified to estimate the imaging PSF which is reparameterized by a Figure 9. The PSFs between the upsampled reference image and the reconstructed HR images of the hilly area.

Further Discussion
The deviation of the PSF between the HR image and the upsampled image is less than the deviation estimated in the PSF estimation stage, i.e., the blur is not removed completely. A visually better image can be generated by contrast manipulation plus sharpening after the SR reconstruction. We take the real Gaofen-4 images of the offshore aquaculture area as an example. The result of enhancing the HR image reconstructed by our SR method (Figure 4c) is shown in Figure 10a. For comparison, the reference LR image is enhanced directly without SR reconstruction, and the result is shown in Figure 10b. The regions denoted by the rectangles in Figure 10a,b are shown in Figure 10c,d. Comparing Figures 4, 5 and 10, we note that the enhanced LR image is visually better than the LR image and the reconstructed HR image, but contains less details. That is to say, the SR algorithm reveals more information contained in the LR images. The PSFs between the HR images and the upsampled reference image of the hilly area are shown in Figure 9. The deviations of the fitted PSFs for IBP, Farsiu's algorithm and our MAP-based SR algorithm are 0.98, 0.82 and 1.0 respectively, which means our method yields the best SR result. Figure 9. The PSFs between the upsampled reference image and the reconstructed HR images of the hilly area.

Further Discussion
The deviation of the PSF between the HR image and the upsampled image is less than the deviation estimated in the PSF estimation stage, i.e., the blur is not removed completely. A visually better image can be generated by contrast manipulation plus sharpening after the SR reconstruction. We take the real Gaofen-4 images of the offshore aquaculture area as an example. The result of enhancing the HR image reconstructed by our SR method (Figure 4c) is shown in Figure 10a. For comparison, the reference LR image is enhanced directly without SR reconstruction, and the result is shown in Figure 10b. The regions denoted by the rectangles in Figure 10a,b are shown in Figure  10c,d. Comparing Figures 4, 5 and 10, we note that the enhanced LR image is visually better than the LR image and the reconstructed HR image, but contains less details. That is to say, the SR algorithm reveals more information contained in the LR images.

Conclusions
In this paper, we propose an SR algorithm to reconstruct an HR image from multiple LR images captured by the Gaofen-4 staring imaging satellite. An appropriate observation model with photometric parameters is built. The registration with photometric adjustment improves the registration accuracy and the quality of the final reconstructed HR image. A MAP-based blind deconvolution algorithm is modified to estimate the imaging PSF which is reparameterized by a

Conclusions
In this paper, we propose an SR algorithm to reconstruct an HR image from multiple LR images captured by the Gaofen-4 staring imaging satellite. An appropriate observation model with photometric parameters is built. The registration with photometric adjustment improves the registration accuracy and the quality of the final reconstructed HR image. A MAP-based blind deconvolution algorithm is modified to estimate the imaging PSF which is reparameterized by a Gaussian function. Finally, the PSF is used to estimate the HR image. Experimental results with synthetic and real images show that the reconstructed HR image reveals more details which cannot be seen in LR images.
Author Contributions: Jieping Xu is the first author and the corresponding author of this paper. His main contributions include the basic idea and writing of this paper. The main contributions of Yonghui Liang and Jin Liu include analyzing the basic idea and checking the experimental results. The main contribution of Zongfu Huang was in providing the real Gaofen-4 images and refining the paper. All authors have read and approved the final manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.