Reference Information Based Remote Sensing Image Reconstruction with Generalized Nonconvex Low-Rank Approximation

: Because of the contradiction between the spatial and temporal resolution of remote sensing images (RSI) and quality loss in the process of acquisition


Introduction
The multispectral sensors of a satellite system sometimes fail or degrade after the system is deployed, and degraded sensors provide blurred and noisy images.For example, one band of the multispectral sensor on board the MODIS satellite degraded after the launch, while the rest of the bands remained unaffected.Damage to images can also be caused by defocusing, atmospheric turbulence and some other factors besides the damage of sensors.Owing to sensor failure and poor observation conditions, remote sensing images are easily subjected to information loss, which hinders our effective analysis of the Earth.Digital imaging processing techniques are commonly used to improve image quality and increase application potential.
In a typical multispectral satellite imaging system, multiple images from different sensors in the same area are available.When one of those images in a multiple image set is degraded, another image in the set can be used as a prior image for reconstruction.However, some image reconstruction approaches have used multiple reference input images.Reconstruction approaches utilize multiple reference input images, which are used as the prior to improve the reconstruction performance in terms of spatial resolution and noise reduction.Li et al. [1] proposed a model utilizing spectral and temporal information as the complementary information for reconstructing the missing information in remote sensing images.Lizhe et al. and Hao et al. [2,3] developed a compressed sensing object function that uses a reference image as a prior.The sparsity constraints in the transform domain come from the target image, and the gradient priors in the spatial domain come from the auxiliary reference image.Peng et al. [4] proposed an algorithm based on the total variation approach using an auxiliary image.These methods have seen the potential for higher reconstruction performance when aided references are present, but the accuracy and quality of reconstruction are lacking.
Since the image prior knowledge plays a critical role in the performance of image reconstruction algorithms, designing effective regularization terms to reflect the image priors is the core of remote sensing images (RSI) reconstruction.To sum up, most of these reconstruction methods are based on some specific prior knowledge of the images, and they use the spectral or spatial information of the remote sensing images to recover the current pixel or to obtain high resolution images.Low-rank-inducing regularization terms have recently received considerable attention in sparse representation.The principle of low-rank approximation is that similar patches are grouped to share a similar underlying structure and form a low-rank matrix appropriately [5,6], and it can be solved by using the efficient singular value decomposition as the optimization tool.Using the low-rank framework, Dong et al. [7] proposed a nonlocal low-rank algorithm, called the spatially-adaptive iterative singular-value thresholding method.Making use of similar information from another relative band, it is easy to extend these nonconvex penalty functions on singular values of a matrix to improve low-rank matrix recovery performance.
In the proposed algorithm, instead of using the 0 norm as the minimization constraint of compressive sensing, the generalized nonconvex low-rank approximation is exploited as the basic unit of sparse representation.The model integrates the wavelet texture feature as structural and textural reference information to establish a novel compressive sensing algorithm for remote sensing image reconstruction (GNLR-RI).Compared to traditional 0 norm-based compressive sensing, the proposed GNLR-RI algorithm makes three contributions: (1) to overcome over-smoothing and texture detail loss of sparse representation, the wavelet coefficient of the reference image is injected as texture reference constraint information; (2) nonlocal similar patches from a reference image are extracted in the low-rank approximation step; (3) the spectral reference information and similarity within bands are jointly used.Therefore, the proposed algorithm can get accurate reconstruction results and improve calculation speed.The solution of the proposed algorithm is derived by the conjugate gradient algorithm, single value threshold (SVT), simultaneously calculating the low-rank matrix of similar image patches and then estimating the reconstructed image.
The remainder of this paper is organized as follows.A brief review is given on the related work of remote sensing image reconstruction with reference images and other generalized reconstruction methods, and auxiliary information is added in Section 2. Section 3 introduces the generalized nonconvex low-rank approximation model, which contains eight surrogate functions.The proposed model of introducing the generalized nonconvex low-rank approximation algorithm and the solution of the proposed model are presented in Section 4. Experiments are illustrated in Section 5. Finally, Section 6 summarizes our conclusions.

RSI Reconstruction with a Reference
Reference images are used to restore missing information of pixel blocks, possibly caused by thick clouds or contrails, invalid detectors, optical dust, etc.For example, toward solving problems with the Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) Band 6 and the Scan Line Corrector (SLC) failure on Landsat ETM+, Li et al. [1] introduced spectral and temporal pixels for reconstructing the missing information in remote sensing images by introducing information from auxiliary images, modulating the injected gray level and regularizing with a synthesis model and an analysis model on the basis of sparse representation.Their experiments showed that the synthesis model is suitable for using both spectral and temporal complementation, while spectral complementation is the best way for the analysis model.
For remote sensing reconstruction or restoration, some work has emerged that uses reference images to improve reconstruction performance.Lizhe et al. [2] proposed a compressed sensing-based model that used an auxiliary image from another sensor or time to reconstruct remote sensing images.Their method assumed that directions of edges and texture are similar in both the reference image and the reconstruction image, and therefore, they built a gradient similarity-constrained cost item to regularize the reconstruction image.Hao et al. [3] proposed a reconstruction model of a remote sensing image using priors of the auxiliary image from another sensor or time.The model employed priors from both the transform domain and spatial domain.Furthermore, the sparsity priors in the transform domain come from the target image, and the gradient priors in the spatial domain come from the reference image.The algorithm is based on the Bregman split method to optimize the hybrid regularization.Peng et al. [8,9] proposed a total variation approach with an auxiliary image to recover a multispectral image.The reference image provides texture and edge similarities to the degraded image by the regulating strength and direction of smoothness in an anisotropic way, and the amount of prior information is adaptively estimated with normalized local mutual information.Similar techniques are implemented in the partial differential equations for denoising of remote sensing images.These methods have the potential for higher reconstruction performance when aided references are present.

Generalized RSI Reconstruction
For ensuring image quality, additional prior information from images is needed.There are generalized remote sensing image reconstruction methods to improve image quality, to increase application potential and to implement digital image processing.The compressed sensing (CS) theory has gained much attention as a fundamental and newly-developed methodology in the information field.The CS theory states that the image that has a sparse representation in a certain domain can be recovered from a reduced set of measurements largely below Nyquist sampling rates [10][11][12].However, Fowler [13] proposed a reconstruction strategy, compressive-projection principal component analysis (CPPCA), which recovers the hyperspectral image dataset using principal component analysis (PCA).Specifically, the CPPCA recovers both the coefficients associated with the PCA transform and an approximation to the PCA transform basis itself.Ly et al. [14] extended the concept of partitioned reconstruction to the spectral dimension of a hyperspectral dataset and incorporated into CPPCA the paradigm of segmented PCA, in which a dataset is segmented into multiple spectral partitions, and PCA is applied in each segment independently.Since multiple predictions drawn for a pixel vector of interest were made from spatially-neighboring pixel vectors within an initial non-predicted reconstruction, Chen et al. [15] proposed a two-phase hypothesis-generation procedure based on the partitioning and merging of spectral bands according to the correlation coefficients between bands to fine-tune the hypotheses.
The transforms that allow the image to have a sparse representation are named sparsifying transforms.One of the most commonly-used sparsifying transforms is finite difference [16], which is more often known as total variation (TV) regularization for image recovery.The popularity of TV regularization lies in its desirable properties, such as convexity, simplicity and its ability to preserve edges.For example, Xiaohua et al. [17] proposed a novel Inverse Synthetic Aperture Radar (ISAR) imaging framework combining a local sparsity constraint and a nonlocal total variation (NLTV).The sparsity is a prior form in which the number of strong scattering points is smaller than that of pixels in the image plane.It plays the role of the classification of the strong scattering points from the clutter background.Xiyan et al. [18] considered the fusion problem as the colorization of each pixel in the panchromatic image and, thus, introduced a term concerning the gradient of the panchromatic image in the function of total variation to preserve edges.Peng et al. [9] indicated that the distributions of pixel intensity of the multimodal image of different CCD sensors differ greatly form each other, and the directions of their edges are very similar.Then, the edge information and these similar structures are used as the important priors or constraints in the total variation image restoration.
On the other hand, similar patches are grouped such that patches in each group share a similar underlying structure to form a low-rank approximation.The low rank induces penalties related to simultaneous sparse coding (SSC), group sparsity and structured sparsity by employing the low-rank characteristic of nonlocal patches in sparse representation.Hongyan et al. [19] introduced a new hyperspectral imagery restoration method based on the low-rank matrix recovery, which suggested that a clean hyperspectral imagery patch can be regarded as a low-rank matrix.Wei et al. [20] proposed a method that integrates the nuclear norm, TV regularization and the 1 norm into a unified framework.The nuclear norm is used to exploit the spectral low-rank property, and the TV regularization is adopted to explore the spatial piecewise smooth structure of the hyperspectral imagery.Usually, nonconvex approaches, like p (0 < p < 1) or the nuclear 1 norm minimization will guarantee a better recovery by directly attacking the 0 minimization problem [21,22].
An important research area is the sparse representation method with dictionary learning (DL), which builds an adaptive basis from particular image instances for sparse approximations [23].Qiegen et al. [16] propose a novel gradient-based dictionary learning method for image recovery that effectively integrates the popular total variation (TV) and dictionary learning technique into the same framework.Specifically, they train dictionaries from the horizontal and vertical gradients of the image and then reconstruct the desired image using the sparse representations of both derivatives.The dictionary learning strategy alleviates the drawback of a fixed basis (finite difference, wavelet, etc.) that a given basis might not be universally optimal for all images [24], but at the cost of non-convexity and non-linearity [25].Most existing DL methods adopt a two-step iterative scheme, where the sparse approximations are found with the dictionary fixed, and the dictionary is subsequently optimized based on the current sparse coefficients [26].
Exact reconstruction of sparse reconstruction can be achieved by nonlinear algorithms, such as orthogonal matching pursuit (OMP) [27], the iterative shrinkage algorithm [28], the Bregman split algorithm [29] and the alternative direction multiplier method [30].Recent advances have shown that replacing the 0 norm with a nonconvex surrogate function can obtain better sparse representation recovery performance.Weisheng et al. [6] proposed a spatially-adaptive iterative singular value thresholding (SAIST) method by revealing the assumption that the basis is orthogonal, and the SSC problem can be rewritten in a low-rank view in some special cases.Shuyuan et al. [31] proposed a multi-task learning and K-SVD-based super resolution image restoration method that learned a redundant dictionary from some example image patches.

Reconstruction Model Based on Reference Images
As in psychological research of the visual cortex, energy distribution in the wavelet coefficient can be identified as the unique features of texture characterization, so we have extra texture features of the reference image and target image.Let α(i, j) denote the wavelet coefficient, F 1_energy denote the 1 norm-based energy, F sd denote the standard deviation and F aad and F entropy denote the average absolute deviation and entropy. (2) where α(i, j) is the wavelet coefficient and ᾱ is the mean value of the wavelet coefficient.Research shows that the combination of these characteristics results in a better performance.Therefore, in this study, the empirical combination with the 1 energy norm F 1_energy , the standard deviation F sd , the average absolute deviation F aad and the entropy F entropy of every wavelet coefficient sub-patch were computed to form the texture feature vector F. As for the reference image, we obtain the texture feature vector F re f in the same way and compute the metric distance to objectively measure the similarity of these two texture feature vectors.The literature shows that not only the texture feature, but also the similarity measurement influence the accuracy of texture feature extraction.Two feature vectors from target images and reference images are as follows: F re f _k = F l1_en α re f _k , F sd α re f _k , F add α re f _k , F en α re f _k (6) where F k denotes the wavelet coefficient sub-blocks of the target image and F re f _k denotes the wavelet coefficient sub-blocks of the reference image.As for the relationship of the kth sub-block of in the reference image and target image, the Canberra distance is adopted to calculate their similarity: In the formulation above, h represents the divergence of these two feature vectors, and the denominator normalizes the divergence.Therefore, the similarity measurement can avoid the scale effects, and the constrained item is formed as follows: where ξ is a small constant, U = M × N/ (m × n) is the number of the sub-blocks of the wavelet coefficients and w k denotes the distance between the wavelet coefficient k-th block of the texture features in the target image and the reference image.More similar texture feature vectors correspond to blocks of wavelet coefficients that are more alike.
The reconstruction object formula of the compressive sensing with the reference information model can be presented as: where y − φψα 2 2 is the 2 data fidelity term and α 1 represents the regularization term denoting prior knowledge.There are many optimization methods to solve the above 1 minimization problem, such as the iterative shrinkage algorithm [28], the Bregman split [29] and the alternative direction multiplier method [30].However, recent advances have shown that replacing the l1 norm with the nonconvex surrogate function can obtain better CS recovery performance.For example, the methods proposed by Fazel [32], Gasso [33], Chartrand [34] and Trzasko [35] have proved that in certain situations, the nonconvex surrogate function is able to recover the sparsity coefficient more efficient.

Generalized Nonconvex Low-Rank Approximation Model
A nonconvex low-rank model for CS recovery exploits the nonlocal structured sparsity via low-rank approximation for image reconstruction.Under the assumption that each exemplar patch of is able to find plenty of similar patches in its neighborhood area, a large number of k-nearest-neighbor searches have been implemented for each exemplar patch in a local window, namely: where T is a pre-defined threshold value and H i denotes the assemblage of positions relating to the patches similar to x i .Under the assumption that these image patches have similar structures, the data-formed data matrix x i has a low-rank property.A data matrix X i = x i 0 , x i 1 , . . ., x i m−1 , X i ∈ C n×m is acquired after the search, which is decomposed into X i = L i + W i , where L i denotes the low-rank matrix and W i denotes the Gaussian noise.Then, the low-rank problem can be solved as: where • 2 F denotes the Frobenius norm and σ 2 ω denotes the variance of additive Gaussian noise.The rank minimization almost is an NP-hard problem; hence, we use the nuclear norm • * (sum of singular values) to replace a series of convex surrogate functions of the rank to obtain an approximated solution.Using the nuclear norm, the rank minimization problem can be efficiently solved by the technique of singular value thresholding (SVT).
In this paper, we consider a smooth, but non-convex surrogate of the rank rather than the nuclear norm.Specifically, according to [32], the rank minimization problem with regard to L i can be approximately solved by minimizing the following function: where • 2 F denotes the Frobenius norm, σ 2 ω denotes the variance of additive Gaussian noise, , and Σ 1/2 is a diagonal matrix whose diagonal elements are the singular values of the matrix L i .By taking a proper parameter λ, Equation ( 12) can be transformed into: For each exemplar image patch, we can approximate the matrix X i with a low-rank matrix L i by solving Equation (13).Many nonconvex surrogate functions have been proposed to achieve a better approximation to the 0 norm [36], including the logarithm function (Log) [37], the p norm (Lp) [38], Geman (Geman) [39], Laplace (Lap) [35] and the exponential type penalty (Etp) [40].In addition, there are several discontinuous functions, such as the smoothly-clipped absolute deviation (Scad) [41], Capped 1 (Cappedl1) [42] and the minimax concave penalty (Mcp) [43].The definitions and super-gradients of these surrogate functions have similar monotonous trends, as displayed in Table 1 and Figure 1.

Nonconvex Surrogate Functions Formulation
Geman

Problem Definition
In multispectral remote sensing applications, the target images' and the reference images' intensity distributions are different, but their edge directions and texture information are often similar.One of the reference images can be used as a prior image in the reconstruction process.This section gives the generalized nonconvex low-rank approximation algorithm for CS recovery exploiting the nonlocal structured sparsity via low-rank approximation for remote sensing image reconstruction.Let x ∈ C N denote the target images, where x is sparsely expressed as x = ψα with a sparse signal α, and y ∈ C M denotes the observed data; the measurement matrix φ ∈ C M×N (M < N) maps x to y.The remote sensing image reconstruction problem is to reconstruct x. Figure 2 shows the relationship between high resolution target images and low resolution observed images.The different low and high resolution images contain different time series or different bands.In some cases, we use one or two different time series or different bands as the reference image and the other as the unknown target image.In this section, a compressed sensing model with reference images integrating low-rank regularization and wavelet textural constraints is proposed and solved with the conjugate gradient algorithm and single value threshold (SVT).Figure 3 shows the framework of the proposed GNLR-RI method.The reconstruction object formula of the generalized nonconvex low-rank approximation with reference information model can be presented as:

High resolution image
where Ri x = R i 0 x, R i 1 x, . . ., R i m−1 x = X i denotes the patch matrix similar to the patch x i .
The variable R denotes the location of the similarity patches; ∑ i R i ψα represents the patch average result.The minimization of the objective constraint function turns into non-constraint optimization.
Target image X band 4

Solving the Proposed Model
As it is difficult to get a closed solution from Equation (15) directly, the optimization process is divided into three steps.
For the first step, the reference image wavelet information is 2 norm regularization; thus, the conjugate gradient algorithm is used to solve it.By setting the derivatives of the objective formula in Equation ( 16) with respect to the sparse coefficient α to zero, we can obtain: (20) and where α = U ∑ j w j α j , we can employ the conjugate gradient algorithm to solve Equation (20).
For the second step, the solution of L i can be obtained by solving the following minimization problem: where Ri x = R i 0 x, R i 1 x, . . ., R i m−1 x = X i denotes the patch matrix similar to the patch x i .By substituting Equation (11) into Equation ( 21), L i can be rewritten as: where σ j is the j-th singular value of G(σ j + ε), which can be solved by using a local minimization method.With the first order Taylor expansion, f (σ) can be approximated as: where σ (k) denotes the value of σ in the k-th iteration.This can be worked out by solving the following equation iteratively, After the constants in Equation ( 24) are ignored, Equation ( 22) is rewritten into: where τ = λ/2η and l + ε denotes the weighted nuclear norm.According to the proximal operator of the weighted nuclear norm, the solution in the (k + 1)th iteration can be obtained as: where U ΣV T denotes the SVT of X i , (x) + = max {x, 0}.For the third step, after the solution of all L i is obtained, the following minimization problem can be solved to reconstruct the wavelet coefficient matrix by using the conjunction gradient algorithm (CG): By setting the derivatives of the objective formula in Equation ( 27) with respect to the sparse coefficient α to zero, we can obtain: Noting that the value of ∑ i RT i Ri represents the number of overlapping patches and ∑ i Ri L i represents the patch average result, we can employ the conjugate gradient algorithm to solve Equation (28).Now, this image reconstruction algorithm is summarized here, is called GNLR-RI and represents the eight nonconvex surrogate functions mentioned above.The detailed description of the proposed method is listed in Algorithm 1.The difference between the NLR-CS-baseline and the GNLR-RI mainly lies in the approximation scheme of the 1 norm and the representation of texture features.Therefore, during the reconstruction process, the low-rank operation strongly encourages the similarity of the patches between the target image and the reference image.

Add the constrained term
to the compressed sensing target objective function.
3. Add the low-rank approximation based on the reference image term.4. Solve the optimization problem Equation ( 20) via the conjunction gradient algorithm.5. Iteration: For (k > K 0 ), do 1) Form a matrix X i making up similar patches of x (k) and set 26) and output α i = α (j) i when j = J.End for End for 6.Solve the optimization problem Equation ( 28) via the conjunction gradient algorithm.End while Output : X = X (K) .

Experiment Results
During the experiments, single-channel and multichannel satellite images from MODIS, Landsat 7, Landsat 8 and Google Earth are used as the simulated data to test the performance of our proposed reconstruction framework.(1) Landsat 7 and Landsat 8 provide PAN images at 15-m spatial resolution and 30-m spatial resolution.(2) MODIS provides images at 500-m resolution.(3) Google Earth provides three-channel color tested images.For the Landsat 7 and Landsat 8 image data, Band 4 with 30-m resolution is the near-infrared band with a better spectrum characteristic in which a body of water has a clear outline, and Band 8 is panchromatic with 15-m resolution; thus, Band 4 is chosen as the target image, and Band 8 is chosen as the reference image, which is down-sampled into 30-m spatial resolution to group similar blocks injected into the reconstructed image to compensate the smoothed areas nonlinearly.For the MODIS data, we use Band 4 as the target image and Band 3 as the reference image to evaluate the proposed algorithm.Furthermore, our proposed algorithm is compared to the well-known reconstruction-based models NLR-CS-baseline [7], orthogonal matching pursuit (OMP) and compressive sampling matching pursuit (CoSaMP) [24] by evaluating the experimental results quantitatively and visually.NLR-CS-baseline uses the standard nuclear norm, which is the l 1 norm minimization.CoSaMp and OMP reconstructs signal by using a very small number of points.Furthermore, our proposed algorithm is compared to the well-known reconstruction-based models NLR-CS-baseline [7], orthogonal matching pursuit (OMP), compressive sampling matching pursuit (CoSaMP) [27], multi-hypothesis predictions combined with block-based compressed sensing with smoothed projected Landweber reconstruction (MH-BCS-SPL) [44], recovery via collaborative sparsity (RCoS) [45] and adaptively-learned sparsifying basis via L0 minimization (ALSB) [46] by evaluating the experimental results quantitatively and visually.NLR-CS-baseline uses the standard nuclear norm, which is the L1 norm minimization.CoSaMp and OMP reconstruct the signal by using a very small number of points.The multi-hypothesis prediction is used to generate a residual in the domain of the compressed-sensing random projections in MH-BCS-SPL.RcoS compels local 2D sparsity and nonlocal 3D sparsity simultaneously in an adaptive hybrid space-transform domain to utilize the intrinsic sparsity of natural images and to confine the CS solution space.ALSB enforces the intrinsic sparsity of natural images substantially by sparsely representing overlapped image patches using the adaptively-learned sparsifying basis in the form of the L0 norm.

Evaluation of the Low Rank Penalty Function and Different Nonconvex Surrogate Functions
First, we justify the necessity of adding the low-rank approximation method.For convenience, the objective function Equation ( 15) recovering the L1 norm without reference information is denoted as L1-WRI.The reconstruction method introduced the reference information constraint with L1 norm Equation ( 9) denoted as L1-RI.The objective function Equation ( 15) that employed the nonconvex surrogate logarithm function is denoted as GNLR-RI-Log.Figure 4 demonstrates the results of these three methods.The proposed algorithm that combines these and introduces low-rank approximation provides better performance.Second, in this section, the reference images of Landsat 7 and Landsat 8 are down-sampled to the same scale of the target images, and the MODIS images are the same scale with 256 × 256, where 256 is the length and width of the input image.Common indices, such as the peak signal to noise ratio (PSNR) and root mean square error (RMSE) are adopted to give a quantitative assessment.Table 2 shows PSNR (dB)/RMSE with different generalized nonconvex surrogate functions.The reconstructed images with different nonconvex functions are displayed in Figure 5.It can be observed that continuous nonconvex surrogate functions, such as Log, Lp, Geman, Lap and Etp, have the same performance as piecewise nonconvex surrogate functions, such as Scad, Cappedl1 and Mcp.In this subsection, we focus on usual nonconvex penalties proposed for recovering sparsity in Figure 1.As all of the penalty functions share common properties, concave and monotonically increasing on [0, ∞), thus their super-gradients are nonnegative and monotonically decreasing.Our proposed general solver is based on this key observation.One of the nonconvex penalties that circumvents Lasso weak points is the SCAD penalty, which has the unbiasedness properties.Among all other penalty functions that lead to sparsity, a popular one is the Lp pseudonorm when 0 < p < 1.The main interest of this penalty resides in its quasi-smooth approximation of the L0 sparsity measure as p tends toward the null value.It can provide sparser solutions than Lasso.For the log penalty, we shift the coefficients by a small quantity to avoid an infinite value when the parameter vanishes.

Performance Comparison for Single-Channel Compressed Sensing
In this subsection, one of the nonconvex surrogate logarithm functions (GNLR-RI-Log) is selected to test the performance of the single-channel experiment.The correlation coefficient (CC) and structural similarity (SSIM) are added to give a quantitative assessment of the reconstruction results.Table 3 shows the performance indexes PSNR (dB), RMSE, CC and SSIM of the Log function with different

Multichannel Reconstruction with References
For the multichannel, Bands 4/5/6 of Landsat 8 data and Bands 1/2/3 of Google Earth data are set as simulated images to be reconstructed in the multichannel reconstruction subsection.Reference images are derived from related Band 8 of Landsat 8 and related gray one of Google Earth.Commonly-used metrics for multichannel images, such as spectral angle mapper (SAM), RASE, relative dimensionless global error (ERGAS) and Q4, are evaluated to check the availability of quantitative remote sensing, since PSNR is not the only index that reflects the reconstruction effect.Note that the ideal result is one for Q4, while it is zero for SAM, RASE and ERGAS.Table 4 shows that GNLR-RI-Log has more advantages with Google Earth data than Landsat 8 data.It can be observed that for Landsat 8 data, our method GNLR-RI-Log works better than NLR-CS-baseline, OMP, CoSaMp and RCoS in terms of PSNR, while it performs worse than MH-BCS-SPL and ALSB.Except ERGAS, for the other three indexes, SAM, RASM and Q4, GNLR-RI-Log obtains comparable and better results than most of the competing methods.On the other hand, for the Google Earth data, GNLR-RI-Log produces more moderate results and performs best among the competing methods in terms of PSNR, SAM, RASM, ERGAS and Q4.

Parameter Evaluation
Similar to the detail-preserving regularity scheme, this subsection evaluates the sensitivity of the proposed method to parameter settings by varying one parameter at a time while keeping the rest fixed at their nominal values.In the reference information constrained reconstruction algorithm, there are two free parameters, λ and η, in the reconstructed object formula Equation (17).Bands 4/8 of Landsat 8 and Bands 2/3 and 3/4 of MODIS data were tested in this subsection.PSNR, RMSE, CC and SSIM values versus the parameters λ and η are plotted in Figure 7.It is obvious that more fine-tuning of the parameters may lead to better results, but the results with the parameter settings are consistently promising.By manually changing the parameters λ and η, some experiments are used to analyzing their variation tendency in Table 5.To make a fair comparison among the competing methods, we have carefully tuned their parameters to achieve the best performance.The parameters of the other competing methods are designed as follow: for OMP and CoSaMp, default parameters (if required as input arguments) are used; for the NLR-CS-baseline algorithm, the main parameters are set as follows: patch size 6 × 6 and similar patches m = 45 are selected for each exemplar patch.For the MH-BCS-SPL method, we use an empirical regularization parameter value λ, and the initial search window is set to w = 1.We have also carefully tuned the parameters of the RCoS and ALSB algorithms for the purpose of hopefully achieving the best possible performance.

Performances with Varied Noise Levels
In this subsection, we discuss the impact of noise on the reconstruction performance of the proposed algorithm with a 0.1 sampling rate, which means observation data in the compressive sensing.One of the nonconvex surrogate logarithm functions (GNLR-RI-Log) is selected to compare to the NLR-CS-baseline, CoSaMP and OMP algorithms.Bands 4/8 of Landsat 8 and Bands 3/4 of MODIS images were tested in the experiments with added Gaussian noise with a mean of zero and variance of (0,4,8,12).As can be seen in Table 6, the proposed GNLR-RI-Log method performs much better than the NLR-CS-baseline method in terms of noise added.Figure 8 shows the reconstruction results for a sample rate of 0.1 and a noise level of σ = 6.It can be observed that our GNLR-RI-Log method can reconstruct the images more precisely than NLR-CS-baseline, CoSaMP and OMP, compensate for the over-smoothness and obtain good group similarities compared to the other algorithms by extracting structural information from reference images.

Computational Complexity
A good remote sensing reconstruction model is expected to be not only effective, but also computationally efficient.The single-channel image is processed on an Intel(R) Core(TM) i7-6700 CPU @ 3.41 GHz with our MATLAB implementation (MATLAB 2015 with 64 bit).The computational cost performances of some representative methods are shown in Table 7, where we use a 256 × 256 single-channel image as the input.As can be seen in Table 7 where we abbreviate iteration to iter., in terms of execution time, reconstruction with GNLR-RI-Log is, as expected, faster than RCoS and ALSB due to low-rank approximation.On the other hand, the execution times of RCoS and ALSB are much slower than GNLR-RI-Log and MH-MS-BCS-SPL due to RCoS and ALSB learning the adaptive sparsifying basis from a fraction of all patches.OMP and CoSamp are the fastest reconstruction methods with their runtimes of 0.106 and 0.031 s per iteration, respectively.In summary, GNLR-RI-Log achieves favorable performance in terms of both reconstruction accuracy and efficiency.

Conclusions
This paper proposes a novel image reconstruction scheme that introduces wavelet coefficients of reference images as reference information based on compressed sensing and generalized nonconvex low-rank approximation.Nonlocal low-rank regularization enables us to exploit the similarity of patches and the nonconvexity of metric rank minimization.In addition, the single value threshold and conjugate gradient algorithms jointly offer a principled and computationally-efficient solution to image reconstruction.This approach has the following characteristics: (1) the wavelet coefficient acted as the texture feature constraint and extracted structural information from reference images to reconstruct remote sensing images; (2) high-resolution images can be achieved via simultaneously using compressed sensing and generalized nonlocal low-rank approximation.Some experiments are made to evaluate the proposed method, and the results show that the proposed method can achieve higher resolution than the state-of-the-art approaches.However, the proposed method has limited performance for the regions and bands where the spatial correlation is not high.In the future, we will further investigate it and use some change detection approaches to improve the performance of our algorithm.

Figure 2 .
Figure 2. Relationship between the high resolution image and the low resolution image.

Table 5 .
PSNR (dB) with different parameters in different images.

Table 7 .
Reconstruction time for a 256 × 256 single-channel image.