Remote Sensing Image of The Landsat 8–9 Compressive Sensing via Non-Local Low-Rank Regularization with the Laplace Function

Utilizing low-rank prior data in compressed sensing (CS) schemes for Landsat 8–9 remote sensing images (RSIs) has recently received widespread attention. Nevertheless, most CS algorithms focus on the sparsity of an RSI and ignore its low-rank (LR) nature. Therefore, this paper proposes a new CS reconstruction algorithm for Landsat 8–9 remote sensing images based on a non-local optimization framework (NLOF) that is combined with non-convex Laplace functions (NCLF) used for the low-rank approximation (LAA). Since the developed algorithm is based on an approximate low-rank model of the Laplace function, it can adaptively assign different weights to different singular values. Moreover, exploiting the structural sparsity (SS) and low-rank (LR) between the image patches enables the restored image to obtain better CS reconstruction results of Landsat 8–9 RSI than the existing models. For the proposed scheme, first, a CS reconstruction model is proposed using the non-local low-rank regularization (NLLRR) and variational framework. Then, the image patch grouping and Laplace function are used as regularization/penalty terms to constrain the CS reconstruction model. Finally, to effectively solve the rank minimization problem, the alternating direction multiplier method (ADMM) is used to solve the model. Extensive numerical experimental results demonstrate that the non-local variational framework (NLVF) combined with the low-rank approximate regularization (LRAR) method of non-convex Laplace function (NCLF) can obtain better reconstruction results than the more advanced image CS reconstruction algorithms. At the same time, the model preserves the details of Landsat 8–9 RSIs and the boundaries of the transition areas.


Introduction
The transmission, reception, and storage of Landsat 8-9 remote sensing images (RSIs) is a critical practical challenge in remote sensing (RS). Therefore, the collection time and processing of massive Landsat 8-9 images have become hot research topics. Compressed sensing (CS) aims to sample/compress the original image using part of the image data (or the corresponding frequency-domain data) and then reconstruct the sampled/compressed data at the terminal by obtaining the reconstructed image [1] close to or beyond the quality of the original data. Generally, CS theory exploits the sparsity and low-rank prior information of the original image to compress the RSI through a data compression method and to recover and reconstruct the compressed data when required.
Nevertheless, the Landsat 8-9 RSIs involve huge real matrices, and this huge data volume (the compressed file size of each group of Landsat 8-9 RSIs exceeds 1GB) imposes a long data transmission time [2], while many practical applications require fast remote sensing data acquisition, limiting the practical application of Landsat 8-9 RSIs. For example, fast data acquisition enables timely assessment of the losses caused by natural disasters, and real-time monitoring of ground objects also depends on the rapid transmission and analysis of RSIs. Therefore, CS technology has prominent practical applications in reconstructing Landsat 8-9 RSIs. At the same time, the compression perception model of an RSI poses a great practical application potential in improving the energy efficiency of imaging sensors [3]. Given that the proportion of random measuring signals acquired by CS in the image is relatively low (typically 10-25%), measuring signal transmission [4] has a significant advantage in time and space storage.
The traditional CS model usually employs a regularization model, compresses the original image using the frequency-domain random sampling (FDRS) method, and then obtains the compressed data [5] with sparsity. Consequently, many image CS reconstruction models are built on practical problem situations. For instance, the traditional CS reconstruction method can be effectively solved using the 11-norm of the sampled signal, i.e., the sparsity of the original signal, and the low-rank function replacement (LRFR) method [6], which effectively solves the CS reconstruction problem. In recent years, the concept of sparsity has been developed into various complex forms, including model-based or Bayesian, non-local sparsity, and structural sparse/group sparsity [7]. Thus, CS data can be reconstructed by exploiting the high correlation between the sparse coefficient, while some experimental studies have demonstrated that the non-convex optimization method (NCOM) based on CS has better image reconstruction results. Indeed, NCOM-based reconstructed images are often higher in visual effect and numerical accuracy than the convex optimization method (COM) [8,9]. To exploit the high correlation between the sparse coefficient and non-convex optimization and find a solution algorithm with relatively high computational efficiency, [10] introduced a method for solving the nuclear norm low-rank approximation model (LRAM) using the SVD algorithm, which yields an easily solved convex optimization problem (COP) [10] by minimizing the sum of all singular values in the image. However, each singular value of the optimization problem in image reconstruction has its practical physical meaning, and therefore, each singular value in any real problem should be treated differently [11]. However, all the singular values have been averaged in the nuclear norm LRAM, limiting the algorithm's ability and flexibility.
Spurred by the findings presented above, we suggest a CS reconstruction model for the Landsat 8-9 RSI based on low-rank approximation (LRA) of non-local associative nonconvex Laplace functions (NL-NCLF), which adaptively assign weights to different singular values [12]. Extensive experimental studies demonstrate that the reconstruction results of NCOMs are more accurate than COMs, but the algorithms are slightly complex and impose a greater computational complexity [13]. In the ongoing big data era, overemphasizing the algorithm's computational time is no longer an important target; therefore, NCOMs have good practical value and prospects. Hence, the core priority is the model to make full use of the prior information of the original signal (such as sparsity and low-rank) in the CS theory of Landsat 8-9 RSIs [14,15].
Generally, the image patches of Landsat 8-9 RSIs present a strong sparsity and have an approximate structural low rank [16]. Therefore, due to the low rank of the structural information in the Landsat 8-9 RSIs, we employ the SVD method to solve the CS model and thus obtain better CS reconstruction results [17]. Moreover, for the image reconstruction, we utilize the regularization model of the low-rank approximate penalty of the non-local combined non-convex Laplace function, which is based on structural self-similarity in the total variation framework [8,9,12]. Our model is a non-convex CS reconstruction model involving a regularization/penalty constraint, low-rank approximation, and image patch grouping [18]. Substituting the low-rank regularization with a non-local and Laplace function enhances the CS reconstruction accuracy compared to common low-rank substitution methods, e.g., nuclear norm, with extensive theoretical analysis proving our algorithm's effectiveness and stability.

Background
In the CS theory of Landsat 8-9 RSIs, random linear measurement/sampling is conducted in the corresponding Fourier transform domain (FTD) of the RSIs x, and then the sampling data y of the Fourier frequency domain (FFD) is obtained. The original Landsat 8-9 RSIs x are recovered through image reconstruction with the sampling data y. In general, the real Landsat 8-9 RSIs are not absolutely sparse but approximately sparse under a certain transformation domain, i.e., they can be transformed into a compressible signal. The sparsity or compressibility of image signals is an important premise and the theoretical basis for CS. Hence, a remote sensing image x can be represented by a linear combination of a set of bases Ψ = {η i }, i = 1, 2, · · · , n, where η i ∈ C N are the sparse bases, i.e., x = Ψα. Specifically, y can be expressed as y = Φx = ΦΨα, where x ∈ C N , y ∈ C M and Φ ∈ C M×N is the measurement/sampling matrix, M < N. Since M < N, matrix Φ is not a full rank, i.e., multiple reconstruction results x ∈ C N can be generated using the same measurement y. Generally, if the measurement matrix Φ satisfies the limited isometry property (RIP) condition [12,16], CS theory can guarantee the perfect reconstruction of signal x by utilizing the sparse (or compressed) signal y. Prior information about the image x must be known to reconstruct a unique perfect image x from the measured data y. The traditional CS method is to recover remote sensing images by using the sparsity of images x, and to satisfy y = Φx, the following constraint optimization problem is employed: where 0 is a pseudo-norm counting the number of non-zero elements in α. However, minimizing the norm 0 is an NP-hard problem. Therefore, we recover the RSI signal from the random measurement/sampling y by using the convex norm l 1 instead of the non-convex norm l 0 and then solve the norm l 1 optimization problem. The optimal model obtained by substituting norm l 1 for the norm l 0 is a convex optimization model, and therefore multiple solution algorithms such as an iterative contraction algorithm (ICA) [19], Bregman splitting algorithm (BSA) [20], and alternating direction multiplier method (ADMM) [21,22] can effectively solve this problem. Recent studies have demonstrated that replacing the norm l 1 with a non-convex norm achieves better CS reconstruction results [16].
By modeling the high correlation between sparse coefficients, the uncertainty of unknown signals can be significantly reduced and afford a more accurate CS reconstruction. The structural sparsity of Landsat 8-9 RSI is particularly important in establishing the CS model for RSIs. Usually, image structure information presents a rich repeatability, so the non-local self-similarity (NLSS) principle of an image structure can be obtained by combining the non-local method [14] and the simultaneous sparse coding (SSC) mechanism [23]. One of the key issues in the CS model is the minimization of the rank function. Because the kernel norm is the minimum convex envelope of the rank function, the rank function in the traditional CS model is often relaxed to the kernel norm (the sum of all singular values of the matrix). However, many works in the literature have proven that using kernel norm to approximate the rank of a matrix has many weaknesses, especially when the matrix has large singular values, the inaccuracy of this approximation is particularly obvious. Therefore, it is particularly necessary to study how to construct a more accurate rank approximation function (RAF) and establish a corresponding low-rank matrix (LRM) recovery model based on it. In recent years, the non-convex function approximation (NNFA) of the LRM's rank function has widely concerned many scholars. A large number of experimental results have proven that these non-convex rank approximation functions are more accurate than the convex function approximation of the kernel norm. The approximate rank function method of the non-convex function (NCF) can not only avoid the NP-hard problem but also provide the optimal solution. The non-convex function replacement (NCFR) can usually provide more scalable solutions. The comparison between the non-convex Laplace replacement function (NCLRF) and the kernel norm replacement function (KNRF) under standard conditions shows that the NCLRF model can better approximate the rank function than the kernel norm model when solving the rank minimization problem. Therefore, it is very meaningful to explore the non-convex regularized low-rank approximation (NCRLRA) model, establish a fast and effective solution algorithm, and apply it to practical problems such as the compressed sensing of a remote sensing image (RSI). This paper suggests a variational framework for CS reconstruction using non-local structural sparsity (NLSS) and non-convex low-rank approximation (NCLRA) [23].

Regularized CS Reconstruction Model Based on Non-Local and Non-Convex Approximate Low-Rank Functions
This section presents a new regularized CS reconstruction model with non-local and an NCLRA, comprising a patch grouping to describe the self-similarity of the images and an NCLRA for low-rank enhancement. Our method assumes that the non-local self-similarity in Landsat 8-9 RSIs is very rich. This assumption implies that for each sample image patcĥ x i ( √ n × √ n) sample patch at position (i), a sufficient number of similar image patches can be found by performing a k-nearest neighbor (KNN) search algorithm in a local window, e.g., 100 × 100, whenx i ∈ C n , namely: where T is a predefined threshold and G i represents the set of locations corresponding to the similar image patch. After the patch is grouped, we obtain the data matrix For each sample image patch x i , each column of X i represents a patch similar to x i (including x i ). Since these patches have similar structures, the data matrix X i formed is low-rank. In practical applications, X i may be interfered with by noise, thus deviating from the ideal low-rank constraint (LRC). Thus, a better representation of the data matrix X i is X i = L i + W i , where L i and W i represent the low-rank matrix (LRM) and the Gaussian noise matrix (GNM), respectively. By solving the optimization problem in Equation (3), the LRM is reconstructed: where 2 F represents the Frobenius norm and σ 2 w is the variance of additive Gaussian noise. Although the low-rank convex substitution method has good theoretical guarantees, the optimization method of non-convex substitution for the rank minimization problem may obtain better recovery results. This paper uses a smooth non-convex function (Laplace function) as the alternative to the low-rank function. Nowadays, several works employ the nuclear norm approximate low-rank function, which provides equal weight to all the singular values in the image patch. However, in many practical situations, the singular values have different physical meanings and should be treated differently [11,12], which is particularly prominent for RSIs. For example, larger singular values represent low-frequency information in Landsat 8-9 RSIs, while smaller values represent high-frequency information and noise. The Laplace function φ(x) = 1 − e −x/ε used in this paper is closer to the pseudo-norm l 0 than the nuclear norm, and thus the sum of the singular values of the Laplace function is closer to the rank function than the nuclear norm. Additionally, the advantage of the Laplace function is that it automatically assigns different weights to each singular value. Based on the above observations, we propose a non-convex low-rank substitution model of the Laplace function, defined as the norm form in Equation (4): where ε denotes a smaller constant value. Note that function X ε is the sum of n singular value functions of the matrix X. So it is smooth and non-convex. Laplace non-convex functions can substitute rank function, which has been proven better when considering information theory [24,25]. On these grounds, the non-convex low-rank approximation optimization model (NC-LRAOM) of Equation (5) is proposed to solve L i .
In practice, the constrained minimization problem can be solved with an unconstrained minimization problem, namely: By selecting the appropriate λ, Equation (6) can be made equivalent to Equation (5). For each sample image patch, an approximate low-rank matrix L i of the matrix X i can be obtained by solving Equation (6), and for each extracted sample image patch the low-rank was enforced on the non-local similar image patch set. Hence, a new non-convex CS reconstruction model is proposed following the proposed low-rank regularization term. The specific CS reconstruction scheme is as follows: whereR i x . = R i 0 x, R i 1 x, . . . , R i m−1 x is a matrix comprising a group of similar patches of each sample image patch x i . The proposed regularized model with non-local combined non-convex approximate low-rank function (NL-NCALF) simultaneously utilizes the group sparsity of similar image patches and the non-convexity of rank minimization, thus being able to obtain better reconstruction results. In the next section, the proposed objective function will be solved effectively using the minimization method of the non-convex function instead of the low-rank function.

Landsat 8-9 Remote Sensing Image CS Reconstruction Algorithm
The proposed CS reconstructed algorithm of the Landsat 8-9 RSI can be solved by minimizing the objective function and the low-rank matrix (LRM) L i of the whole image x. For the initial estimate of the unknown image x, we first extract the sample image patch x i at each pixel i in each direction and assign a similar set to each image patch x i , as described in Section 3. Then we solve the following minimization problem for each image patch as follows: Equation (8) is solved using Theorem 1: Theorem 1. Given Z ∈ R m 1 ×m 2 , the minimum value of Equation (9): is given by a weighted singular value threshold, i.e.,: ε denotes the gradient of φ at position σ i , and σ i is the ith singular value of X.
Proof. Equation (4) shows that the function X ε is the sum of n singular value functions of matrix X, then Equation (9) can be expressed in the form of Equation (11): where X, Z ∈ C m 1 ×m 2 , and the optimal solution of Equation (11) can be obtained using the general weighted singular value threshold [26,27], namely: Consequently, the conclusion of Theorem 1 can be confirmed.
In the case of a real matrix, the proposed CS reconstruction model for the Landsat 8-9 RSI is non-convex, so it is not expected to find its global minimum. However, the weighted singular threshold algorithm can obtain the minimum value (possibly the local minimum value) of Equation (9). Note that even though the weighted threshold method is only a local minimum, the value of the objective function still decreases. In the experiment, we set w (0) = [1, 1, . . . , 1] T .
According to [28], weighting l 1 -norm performs much better than l 1 -norm in approximating l 0 -norm and generally produces better image CS perception reconstruction results. Similarly, the experimental results in the next section demonstrate that the Laplace function produces CS reconstruction results better than the nuclear norm. By solving the following minimization problem to obtain L i , we can further reconstruct the whole Landsat 8-9 RSI: when the measurement matrix Φ is the Fourier transform matrix (the Fourier transform is one of the important transforms in the field of Landsat 8-9 RSI processing), Equation (12) can be quickly solved using the alternating direction multiplier method (ADMM) [29]. In this case, first, the augmented Lagrangian form of Equation (12) can be expressed as in Equation (13): where z ∈ C N is an auxiliary variable, µ ∈ C N is a Lagrange multiplier, and β is a positive constant. The advantage of ADMM is that Equation (13) can be divided into two subproblems, where both can find their closed solutions. Applying ADMM to Equation (13), the iterative form of the following Equation (14) can be obtained: In Equation (14), ρ > 1 is a constant. By determining x (l) , µ (l) , β (l) , and z (l+1) , a closed solution can be obtained: Note that ∑ iR T iR i is a diagonal matrix. Therefore, Equation (15) can be easily calculated. The sub-problems of x and y can be calculated using Equation (16): where Φ is the Fourier transform matrix (FTM), Φ = DF and where D and F represent the down-sampling matrix and Fourier transform matrix (FTM), respectively. In other words, the sub-problems of x and y can be solved from the image space to the Fourier frequency domain space by Equation (16). Equation (17) can be obtained by replacing Φ with DF in Equation (16) and executing Fourier transform on both sides of the equation.
So, we simplify Equation (17) and calculate x by taking the inverse Fourier transform, namely: By updating x and z according to Equation (14), µ and β can be easily calculated.
After the unknown image x is calculated, the low-rank matrix L i is updated using Equation (8), from which we obtain the updated L i that is used to recalculate the estimate of image x. This process is repeated until the algorithm meets the convergence condition.
The proposed CS reconstruction model for the Landsat 8-9 RSIs is non-convex. Therefore, the hot start method is used to preprocess the CS reconstruction images. When the reconstruction results reach a certain accuracy, the above ADMM method is used to solve the proposed model and obtain higher accuracy. To save computing time, image patch grouping is not updated after each iteration but at the end of an outer cycle.

Data Sources
This study employs Landsat 8-9 RSIs for all simulation experiments obtained from the US Geological Survey Center for Earth Resource Observation and Science (EROS). Landsat 8 is the eighth satellite of the US Landsat Missions (Landsat), successfully launched on an Atlas-V rocket from Vandenberg Air Force Base, California, on 11 February 2013, originally known as the Landsat Data Continuity Mission (LDCM). Landsat 8 carries an operational land imager (OLI) and a thermal infrared sensor (TIRS). The OLI measures the spectrum's visible, near-infrared, and shortwave infrared portions (VNIR, NIR, and SWIR). The TIRS measures the temperature of land surfaces in two thermal bands with a new technology that applies quantum physics to detect heat. The OLI Land Imager includes nine spectral bands with 30-m multi-spectral spatial resolutions, which include a 15-m panchromatic band. This study's first set of satellite images was the Landsat-8 LITP product of the blue band (30-m) with a wavelength range of 450-515 nm. Landsat 9 is the ninth satellite of the US Landsat Missions (Landsat), successfully launched from the Vandenburg Space Force Base, California, on 27 September 2021. Landsat 9 Carries the second-generation land imager (Operational Land Imager 2, OLI-2) built by Ball Aerospace & Technologies and the second-generation thermal infrared sensor (Thermal Infrared Sensor 2, TIRS-2) built by NASA Goddard Space Flight Center. The OLI-2 captures images of the Earth's surface in visible, near-infrared, and shortwave-infrared bands, increasing the radiative measurement accuracy from 12 to 14 bits with Landsat 8 and slightly improving the overall signal-to-noise ratio. The OLI-2 land imager includes nine spectral bands with a 30-meter spatial resolution, which includes a 15-meter panchromatic band. The TIRS-2 measures the thermal infrared radiation or heat of the Earth's surface in two bands that perform better than the thermal band of the Landsat 8. This study's second set of satellite images is the Landsat 9 L2SP product of the red band (30-meter) with a wavelength range of 640-670 nm. Table 1 summarizes the locations and dates of the observed images. The Landsat 8 RSIs were ordered in four different geographical regions: Antarctica; Selkirk in Manitoba, Canada; Flathead Lake in Montana, USA; Lincoln in Washington, USA. Antarctica is mainly composed of ice sheets, glaciers, and snow, while Selkirk in Manitoba, Canada, is a grassland area with relatively flat terrain, a Precambrian shield, and its northern end is the permafrost layer. Flathead Lake in Montana, USA, is in the state's northeast and is known for its rich rocks and plains. It includes rich land features, such as water, mountain, forest, and vegetation. Lincoln has abundant landforms, abundant rainfall, and large desert areas in the east.
The Landsat 9 remote sensing images (RSIs) [30] were ordered in five different geographical regions: Lake Abitibi in Ontario, Canada; Yeosu, South Korea; Shanghai, China; Huangshi, Hubei, China; Erenhot, Inner Mongolia, China. The Lake Abitibi waters range from western Lake Forest to the St. Lawrence River in eastern Cornwall. Lake Abitibi stretches between Ontario and Quebec and includes many land features, such as lakes, forests, cities, and farmland. Yeosu is located in the Yeosu Peninsula, the southernmost tip of the Korean Peninsula, belonging to South Jeolla Province, and it includes the ocean, coastline, and city. Shanghai is located on the west coast of the Pacific Ocean, where the Yangtze River and Huangpu River converge. Its land features include rivers, oceans, cities, islands, and peninsulas. Huangshi City is located on the south bank of the middle reaches of the Yangtze River, facing the Yangtze River in the northeast and Huanggang City across the river, and it is rich in minerals. In the northwest and middle of the city, lakes such as Baoan, Sanshan, and Daye form a plain lake marshland, while the rest are low mountains and hills. Erenhot is flat and gently inclined from southwest to northeast. The land features include forests, deserts, grasslands, villages, and vegetation.
The objective of selecting Landsat 8-9 RSIs of different regions is to ensure that the newly-developed model performs very well for diverse land cover types.

Experimental Results and Analysis
To demonstrate the effectiveness of the newly developed method, extensive simulations were conducted using Landsat 8-9 RSIs. Since Landsat 8-9 RSIs (30-meter resolution) are too large, we adjusted the images by cropping, scaling, recombining, and synthesizing. After processing, the simulated Landsat 8-9 RSI had a 1000 × 1000 resolution used to evaluate the proposed CS model. Then, the CS model was applied to reconstruct the compressed sampling data and obtain the restored image, which was compared against the original image to assess the model. All simulations were conducted on an Intel (R) Core (TM) i9-10980XE CPU @ 3.00 GHz and with 128 GB RAM. The synthetic RSIs for reconstruction are illustrated in Figure 1:  Figure 1 illustrates a series of synthetic images obtained from the original Landsat 8-9 images, which were used to evaluate our algorithm after cutting and resizing. The above images were used in the newly developed CS model for reconstruction and comparative evaluation against other CS models, namely, NL-Laplace-CS, NL-SRF-CS [16], KCS-GSR [31], and NLDR-CS [32], and revealed the advantages and disadvantages of each model. The Landsat 8-9 RSIs were sampled using different frequency-domain sampling ratios (10%, 15%, 20%, and 25%) and were reconstructed using the above model. The corresponding reconstruction results are presented in Figures 2-10. Figure 2 reveals that the NL-Laplace-CS model performs better than NL-SRF-CS, KCS-GSR, and NLDR-CS considering reconstruction. From the first line (10%), the fourth line (15%), the seventh line (20%), and the tenth line (25%), it is evident that the images reconstructed by the NL-Laplace-CS model have no obvious edge blur, but the competitor method reconstructs images with a visible edge blur. The reconstructed image using the NL-SRF-CS model has information loss in surface features and heterogeneous regions of textures. The images reconstructed with the NL-Laplace-CS and KCS-GSR methods present fewer zigzag effects than the other models. The overall reconstruction result of the NLDR-CS model is better than the NL-Laplace-CS and KCS-GSR models but still slightly inferior to the NL-Laplace-CS model. The NL-Laplace-CS model has inherent noise suppression and, therefore, can suppress areas of high-frequency changes. However, from the different images in Figure 2e red magnified diagram at the lower left corner of each image), we observe that the NL-Laplace-CS model retains the most details. Figure 3 illustrates the CS reconstruction results of the Landsat 8 RSI covering Selkirk, Manitoba, Canada, using the same model presented in Figure 2 Compared with the competitor models, the NL-SRF-CS and KCS-GSR models present large differences compared to the original images. Specifically, the images reconstructed with the NL-SRF-CS method have significant edge blurring and lost more details in areas with more edges. The images reconstructed with the KCS-GSR model have a sharp zigzag effect. Moreover, the results of images reconstructed by NLDR-CS show a greater data loss of surface features and edge areas than those reconstructed images using the NL-Laplace-CS model. The overall visual effect of the reconstructed image from the NL-Laplace-CS model has the highest fidelity relative to the original images. Figure 4 depicts the CS reconstruction of a Landsat 8 RSI of Lake Flathead, Montana, USA, using the same model as in Figures 2 and 3. The difference between the reconstructed image of the KCS-GSR model and the original image is the largest, followed by the NL-SRF-CS model and NLDR-CS model. However, the proposed NL-Laplace-CS model has the smallest difference from the original image. The reconstructed images of the NL-SRF-CS model have much edge blurring and loss of detail in the edge area, and the KCS-GSR reconstructed image also lost more details than the original image. The reconstructed image of the NLDR-CS model affords fine surface structure and edge information, but there is a minor blurring effect in the edge area. In the image reconstructed by the NL-Laplace-CS model, the partially magnified area (lower left red box diagram) reveals a good visual effect, while the details of edges and non-smooth areas are also well preserved. Figure 5 illustrates the resulting images of CS reconstruction of Landsat 8 RSIs from Lincoln, Washington, USA, using the model of Figure 2. The difference images show that the reconstructed result images of the KCS-GSR model and the NLDR-CS model differ the most from the original image, and the proposed NL-Laplace-CS model differs the least. For KCS-GSR and NLDR-CS models, the reconstructed image's local magnified area (red block diagram in the lower left corner) has a more obvious edge blur. Furthermore, the reconstructed images of the NL-SRF-CS model reveal that the reconstruction results are significantly serrated, and there is also a small fuzzy effect in the uneven areas. For the NL-Laplace-CS model, the reconstructed images afford a good visualization, the suppression of details is invisible, and the boundaries of each object region in the image are well preserved.     Table 2 summarizes the number of pixels within a range (±260) of around 0 for each image. The more pixels within a pixel value range around 0, the better the reconstruction effect. Table 2   To further evaluate our model's performance, we employ three quantitative image quality indicators (PQI): root mean square error (RMSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM) [33]. A lower RMSE value, higher PSNR value, or higher SSIM value indicates a better model performance. Table 3  Moreover, the average error generated by the NL-Laplace-CS model in all geographic regions is the smallest among all models for all sampling ratios. For Selkirk in Manitoba, Canada (Figure 3), all models have similar PSNR values, but the PSNR values of NL-SRF-CS and KCS-GSR models are relatively smaller than the other models. Although the NLDR-CS model has high PSNR values, it changes relatively slowly with the sampling rate increasing. In any case, the PSNR of the NL-Laplace-CS model is the maximum. Furthermore, the PSNR results (see Table 3) indicate that the NL-Laplace-CS model performs well, and its SSIM value is greater than the competitor models. Although the SSIM value increases with the image sampling rate, the rate at which the new model increases is much larger than the other models. These PQI results indicate that the newly developed NL-Laplace-CS model has lower RMSE values, higher PSNR and SSIM values, better visual effects, and edge preservation. The second group of satellite images used in our simulations involves the Landsat 9 L2SP product in the red band (30-meter) with a wavelength range of 640-670 nm. Similar to the Landsat 8 synthetic simulation images, the Landsat 9 image was trimmed and scaled to obtain the simulation Landsat 9 data presented in the second row in Figure 1. Finally, we evaluated that dataset with the NL-SRF-CS, KCS-GSR, NLDR-CS, and NL-Laplace-CS models and analyze the performance of each algorithm through the simulation results. Figure 6 illustrates the CS reconstruction images of the Landsat 9 RSIs from Lake Abitibi, Ontario, Canada, using the model presented in Figure 2. In the algorithm simulation experiment in this region, various CS reconstruction models can reconstruct the resulting images well. However, there are still some differences, as Figure 6 highlights the fact that the NL-SRF-CS model has a good reconstruction effect in the relatively smooth area, but the reconstruction effect in the edge area and the detailed area is unsatisfactory. The reconstructed resulting images of the KCS-GSR model are quite different from the original images; namely, there is a certain gap between the overall reconstruction results of the KCS-GSR model and those of the NLDR-CS model and NL-Laplace-CS model. For NL-Laplace-CS and NLDR-CS models, the reconstructed image's local magnified area (red block diagram in the lower left corner) has more obvious detail preservation. At the same time, the reconstruction results of the proposed NL-Laplace-CS model have the least difference from the original images. Figure 7 depicts the resulting images of the CS reconstruction of the Landsat 9 RSI of Yeosu, South Korea, using the model presented in Figure 2. For this area, the reconstructed images of all CS reconstruction models present unsatisfactory results, mainly due to the geographic details in this region. Figure 7 reveals that the KCS-GSR model's reconstructed images differ the most from the original ones, followed by the NL-SRF-CS model. The NLDR-CS model shows better reconstruction results but is still imperfect in the marginal and detailed areas. The reconstruction results of the proposed NL-Laplace-CS model present minor differences from the original images and thus affords the best visual effect. Figures 8-10 demonstrate each competitor algorithm's advantages and disadvantages based on each image's corresponding reconstruction result. However, due to paper length limitations, we neglect a detailed analysis. Additionally, Table 4 summarizes the number of pixels within a range (±260) of around 0 for different Landsat 9 images. The more pixel values within the range around 0, the better the reconstruction effect. Table 4 reveals that most pixels of the NL-Laplace-CS model-based reconstructed images are concentrated within the range [−260, 260], as illustrated in Figures 6-10. For 16-bit images (0 to 65,535), the range of (0 to 260) is very small. Table 4 highlights the fact that the NL-Laplace-CS model has the largest percentage of pixels (percentage of the whole image) in a range of around 0, indicating that its performance is best among the models evaluated.
To further evaluate the proposed model's performance, we calculate the RMSE, PSNR, and SSIM values on the CS reconstructed Landsat 9 images. Table 5 summarizes the information presented in Figures 6-10 and highlights the fact that the NL-Laplace-CS model has the lowest RMSE value and the highest PSNR and SSIM values. The resulting images used for the CS reconstruction based on the NL-Laplace-CS model had the best visual effect, namely, the most preserved edge information and image details.

Discussion
The NL-Laplace-CS regularization method relies on the appropriate adjustment of the data fidelity and weight parameters in the regularization terms. With the optimal choice of the weight parameters, this method yields the best reconstruction results. The experiments use the Landsat 8-9 RSI for CS simulation reconstruction, and we compared the reconstruction results obtained from the Landsat 8-9 RSI simulation with the NL-Laplace-CS model against the NL-SRF-CS, KCS-GSR, and NLDR-CS image reconstruction models. We also compared these reconstruction results with the original Landsat 8-9 RSI. The results highlighted that the new NL-Laplace-CS had lower RMSE values, higher PSNR and SSIM values, and better visual effects than the competitor models. Errors mostly occur between geological structures and land cover classes in the transition zone. The resulting images of the CS reconstruction of the newly developed NL-Laplace-CS model reveal that, even in such a transition zone, the reconstructed pixel values still have a high fidelity relative to the original Landsat 8-9 RSIs. This is important because most remote sensing image applications require high fidelity in the transition zone to clarify the geological surface features with complex geographical compositions.

Conclusions
This work developed a Landsat 8-9 RSI CS reconstruction model based on a non-local framework combined with low-rank regularization approximation of non-convex Laplace functions, which was solved using the ADMM algorithm. The proposed NL-Laplace-CS model was challenged against advanced CS image reconstruction models such as NL-SRF-CS, KCS-GSR, and NLDR-CS in a simulation reconstruction for Landsat 8-9 RSIs. The simulation results revealed that the proposed NL-Laplace-CS model effectively utilized similar image patches of sparse groups and established a low-rank regularized approximate minimization model of non-convex Laplace functions, suggesting a very effective method for CS reconstruction of Landsat 8-9 RSIs. Specifically, the proposed NL-Laplace-CS model has lower RMSE values, higher PSNR and SSIM values, and outperforms state-of-the-art CS image reconstruction models such as NL-SRF-CS, KCS-GSR, and NLDR-CS.  Institutional Review Board Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.