Abstract

Incomplete remote sensing image observation data require proper processing to restore the original data. Nevertheless, the lack of information, due to noise and dead pixels, along with the quality of the acquired images that may be severely degraded, limits the accuracy of later processing. Therefore, this paper proposes a remote sensing image restoration scheme that employs the low-rank structure of the Bayesian probability tensor latent space and a new tensor completion method, which has strong robustness for the model selection process of the following remote sensing image data analysis. Considering the data space, we incorporate an adaptive nuclear norm regularization into the processing of latent factors that performs singular value decomposition on a more accurate scale. The proposed model adopts two algorithms: Augmented Lagrange Multiplier (ALM) and Alternating Least Squares (ALS). Several synthetic and real-data-based experiments illustrate the inherent ability of our method in restoring original information and preventing overfitting problems, even when a vast number of items are missing. In addition, the combination of adaptive nuclear norm regularization and tensor decomposition is stable, settling a high sensitivity of rank selection and low computational efficiency in traditional tensor reconstruction methods.

1. Introduction

Remote sensing instruments have aroused widespread concern in the research field, as these are highly appealing for observing, learning, and understanding Earth. Nevertheless, severe weather conditions and sensor failures cause the remote sensing images to be heavily disturbed during the acquisition process, significantly reducing their visual quality and affecting the accuracy of the applications where these images are applied [1, 2]. Consequently, extracting raw data from remote sensing images containing dead pixels, clouds, and noise is a significant study subject [35]. However, to exploit efficiently raw remote sensing imagery, it is necessary to forecast and reconstruct missing information [6].

In recent years, advanced technologies such as wavelet compression and anisotropic diffusion have been applied in remote sensing image restoration, yet they neglect the relativity between spatial pixels resulting in poor postprocessing results. Literature offers numerous relevant studies, such as proposing an effective nonlocal method [7] that analyzes the structural similarity through the gray correlation of coefficients and then setting the similarity weight function. Li et al. used a multitemporal dictionary learning algorithm to recover remote sensing images polluted by large thick clouds and the consequent cloud shadows and verified the effectiveness of the two methods [8]. Lorenzi et al. propose three alternatives for reconstructing missing regions in images acquired by very high spatial resolution (VHR) sensors based on the inpainting method [9]. Zhuang and Bioucas–Dias develop a very sparse and compact Hyper-Spectral Imagery (HSI) representation using HSI's low-rank self-similar feature [3]. Xue et al. proposed a new HSI noise reduction method based on CP decomposition modeling with tensor nonlocal LR and automatic rank determination. The method considers the a priori knowledge of GCS and NSS in space and combines the rank estimation of CP decomposition. The results show that the method can effectively reduce the complex noise and maintain the fine weave structure of HSI [10]. Zeng et al. proposed an HSI recovery model based on local tensor kernel parametric minimization and global L1-2SSTV regularization for the HSI recovery problem. The results of extensive experimental simulations show that the model has a good quality evaluation index and visual effect of model recovery results [11]. Dian and Li proposed a subspace-based low-tensor multirank regularization method to improve the spatial resolution of HSI. Finally, the superiority of the method is demonstrated by experimental results on simulated and real datasets [12]. Yang et al. proposed a variational depth HSI-MSI fusion model, and then expanded the fusion model into a DL network, in which the degenerate model and the prior data can be represented by the DL network, and validated the effectiveness of the method from experiments on several HSI datasets [13]. The purpose of both methods is to eliminate one or two types of noise. However, the acquired remote sensing images are generally mixed with various noise types, i.e., edge destruction, dashed noise, impulse noise, and Gaussian noise, and therefore a solution restoring the image and enhance its quality regardless of the noise type is mandatory.

Low-order tensors are an effective means to describe the inherent geometric characteristics of remote sensing images. Some methods involve low-rank matrix modeling, e.g., Zhao and Yang exploit the low rank of high-order matrices and establish a low-rank matrix approximation based on small blocks [14]. Hao and Su suggest noise adjustment and propose an iterative low-rank matrix approximation denoising method, which uses the nonlocal self-similarity of the image. This method adopts a new patch-based model to restore the low-rank tensor applied to denoise multiframe images [15]. Image restoration is essentially the process of estimating the original image, and its purpose is to obtain the optimum estimate of the original image under specific objective standards. The prior knowledge of the image degradation process is related to the effect of image restoration. For reconstructing remote sensing images, low-rank tensor completion is one of the most effective methods to integrate prior knowledge [15, 16]. For instance, both the Tucker and CandeComp/Parafac (CP) models rely on prior knowledge of the model’s original multidimensional data [17]. The CP decomposition can be solved by Alternating Least Squares (ALS), that is, updating a single variable by other fixed variables at each round of iteration. However, this model is nonconvex, subject to local minima, and requires good initialization to get better results [18]. For fixed n-rank, the uniqueness of Tucker decomposition cannot be guaranteed, and some constraints are generally added, such as orthogonality constraints on the factor units obtained from the decomposition [19, 20]. Therefore, we need a new tensor completion method to enhance the robustness of the model selection process for subsequent remote sensing image data analysis.

Low-rank tensors utilize rank minimization to estimate low-rank approximations of incomplete data [2123]. However, low-order tensors may cause overfitting problems due to neglecting critical information of the original data. The Bayesian probability inference is used to overcome this, as the Bayesian probability model requires fewer observation data to form inferences. It can overcome extremely scarce conditions and can be effectively combined with tensor decomposition [24, 25]. Consequently, the matrix/tensor decomposition probability model has received widespread attention in collaborative filtering and matrix/tensor complement. Hence, [26] proposes the probability matrix factorization combined with the Markow Chain Monte Carlo (MCMC) inference to perform complete Bayesian processing. Babacan [27] uses variational Bayesian inference to perform complete Bayesian processing on the tensor structure, while [28] further expands the nonparametric variables and robust variables. Chu and Ghahramani [29] represent the probability framework of tensor decomposition, with other variants involving the extension of the nonparametric Bayesian model [30] and the exponential family model [31].

This article introduces a complete Bayesian probability tensor decomposition model that is appropriate for the TR decomposition framework. Precisely, based on the predicted distribution, incomplete tensor, and the noise of missing items, a potential multilinear factor is calculated to implicitly and automatically determine the rank of the real tensor. To this end, a sparse-induced hierarchical prior based on multiple factor matrices is proposed, and the hyperparameters of each factor matrix are related to each potential dimension so that the elements in the factor matrix are limited and minimal. All model parameters are considered variables that can be assumed, including noise accuracy and corresponding priors. Aiming at the shortcomings of inherent high sensitivity and low efficiency of rank selection of traditional tensors, we design a combination of adaptive tensor nuclear normalization methods. The model comprises Bayesian probability tensor decomposition and an adaptive kernel model and affords good performance, robustness, and high computational efficiency. The Alternating Least Square (ALS) and Augmented Lagrange Multiplier (ALM) algorithms are adopted to optimize this model.

This paper primarily studies the problem of remote sensing image reconstruction under the influence of mixed noise. The main contributions of this work are summarized as follows:(i)A Bayesian probability tensor decomposition method is used to fully exploit the global space information heavily affected by complex noise. This strategy affords a low computational burden, and the tensor decomposition model is constructed using actual observation images.(ii)The designed adaptive nuclear norm regularizes the tensor’s decomposition rank, affording our algorithm to obtain a stable solution. The results manifest that the regularization enhances the remote sensing image information and better reserves the specific information that remote sensing images capture.(iii)The Augmented Lagrange Multiplier (ALM) and Alternating Least Square (ALS) algorithms are employed to solve the adaptive nuclear norm regularized Bayesian probability tensor decomposition model. A considerable amount of simulations and actual data experiments demonstrate that the algorithm is more effective than current commonly used methods, both in visual comparison and quantitative evaluation.

This essay is structured as follows: Section 2 introduces the Bayesian Gaussian TR decomposition and the remote sensing image reconstruction problem formulation, which will introduce our model. In Section 3, the adaptive nuclear norm regularized Bayesian probability tensor decomposition model is introduced, and the ALS and ALM algorithms are adopted to solve the proposed model. Extensive experimental simulations of both methods are implemented in Section 4; while finally, Section 5 concludes this work.

2. Bayesian Gaussian TR Decomposition and Problem Formulation

2.1. Notations

In this paper, the notations of Kolda and Bader [32] are utilized, i.e., tensors of order are represented by , matrices by , vectors by , scalars . represents the element of tensor at position . Furthermore, two kinds of tensor expansion processing types are adopted, where the first mode-n expansion [31] of tensor is represented by , and the second one, which is usually adopted in Symmetric tensor decomposition processing, is represented by .

2.2. Bayesian Gaussian TR Decomposition

Tensor decomposition is widely used to reduce the storage cost and operational complexity of high-dimensional data and involves a linear scale transformation of the input dimension instead of an exponential scale transformation. refers to a higher-order tensor, where is the size of . The TR decomposition of can be considered as follows:where denotes the -th element of the tensor and is the -th lateral slice matrix of the latent tensor .

Regarding the Bayesian probability model, it is assumed that the noise contained in the observed entry obeys a Gaussian distribution:where represents the Gaussian distribution and is the accuracy.

To better fit the tensor data, a flexible prior part is placed on and , and it is assumed that this prior distribution applies to all factor matrices and is multivariate Gaussians:

To boost the model’s robustness and efficiency, the conjugate Gaussian-Wishart priors for hyperparameters and are placed [33].where is the degree of freedom of the Wishart distribution and is given by the following equation:

The entire graphical structure characterizing the Bayesian probabilistic TR decomposition data generation process is illustrated in Figure 1. Given that the model has a ringlike graphical appearance, it is intuitively called tensor ring decomposition. In this graph model, is the observation data, and are the parameters, and and are the superparameters. For further information on setting these superparameters when initializing the sampling algorithm, the reader is referred to [23].

2.3. Problem Formulation

In many practical situations, the observed remote sensing image will be contaminated by various noise types [25, 34, 35]. Hence, a remote sensing image y can be represented by a linear regression model:where and have the same size representing the observed vertical image and the mixed noise items, respectively. The purpose of rereflecting the image is to predict the complete reflection of image X based on the observed incomplete data Y.

In this work, the Gaussian noise term and the sparse noise term are divided from noise E. Considering that the sparse noise term includes stripes and impulsive noise:

Therefore, the Frobenius norm and norm can be used naturally to simulate these two noise terms and , respectively.

3. Remote Sensing Image Restoration

3.1. Proposed Model

The purpose of tensor decomposition is to search for latent factors in the tensor data, which can be regarded as latent features of tensor data. The observation object is represented by a low-rank structure that is generated through tensor decomposition. According to the definition of tensor decomposition, an optimization model is drawn up:where represents all observed entries, and a set of observation indices is . denotes the approximate tensor produced by the latent factor and a potential factor generated by the Bayesian probability TR tensor decomposition.

In particular, the Bayesian probability tensor factorization model adopts a hierarchical probability model to represent the TR factorization. Specifically, several potential factors are introduced into the space induced prior, the complete Bayesian method is adopted, and the appropriate priority is introduced to determine the rank. A remote sensing image can be combined with several noise types, like fringe, Gaussian noise, dashed line, and pulse noise, to simulate real-world scenarios. The bayesian probability tensor decomposition model obtains the suboptimal rank, and consequently, the rank minimization tensor complement model is introduced:where is the low-rank tensor after recovery. Therefore, the model can discover the low-rank structure of the data and approximate the recovered tensor. Additionally, an adaptive tensor nuclear norm is designed as follows:where is the weighted nuclear norm of and is a nonnegative weight related to the Gaussian noise variance in each set of observations. Depending on the importance of the image center point or field, different positions are given different weights. Considering the weight of the edge position of the observation data as one, while the internal position increases proportionally, the center position has the most prominent weight to ensure that each weight is an integer and reduces the amount of calculation. Embedding this in the model of rank minimization tensor completion,

3.2. Optimization

To optimize the proposed solution, we employ the ALM and ALS algorithms shortly introduced as follows.

3.2.1. ALM

Equation (11) is rewritten as an equivalent minimization problem by introducing some auxiliary variables:where is the weighted N-dimensional difference operator. The problem of equation (12) is transformed into the minimized augmented Lagrangian function of equation (13), under the constraints , where is the penalty parameter, and are the Lagrange multipliers. Consequently, the augmented Lagrangian function (equation (13)) can perform alternative optimizations to one variable while modifying other variables. Variables involved in the model of equation (8) can be updated in the k-th iteration as follows:

From the augmented Lagrangian function (13), all terms containing are extracted as follows:

Equation (14) can be transformed into the following equal problem:

Parameters and can be calculated using the classical Bayesian tensor algorithm. Then,

Considering also the augmented Lagrangian function (equation (13)), can be updated towhere represents the adjoint operator of . The block-circulant structure of the matrix corresponding to the operator can be diagonalized by the FFT matrix, and thuswhere fftn and ittfn represent the fast Fourier transform and the inverse fast Fourier transform and is the elementwise square. The division is also implemented elementwise.

By extracting from equation (13) all terms containing ,

Moreover, by introducing a threshold contraction process,where represents the threshold shrinkage processing.

Similarly, by updating ,

Using the threshold shrinkage processing introduced above, the following formula is obtained:

Then, by updating through extracting from equation (8) all items related to ,

Finally, the following formulas are used to update the multipliers based on the ALM algorithm:

3.2.2. ALS

ALS is adopted to settle the model represented by equation (13) and includes two steps: kernel tensor update and adaptive nuclear norm regularization reconstruction.

Initially, the core tensors are updated sequentially. For n = 1, …, N, the optimization iswhere is the elementwise multiplication. We adopt , , , and to denote the matrices , , , and , respectively, where the core can be updated on a row-by-row basis. Let be the row index of , then the update of iswhere is a vector extracted from the i-th row and the observed column provided by . To resolve equation (26), the least square method is adopted, and the optimal solution is given by the following equation:

Subsequently, once all core tensors are updated, the adaptive nuclear norm regularization is adopted to reconstruct tensor .

4. Experiments and Discussion

The validity of the proposed method is demonstrated through real-data and simulated experiments.

4.1. Simulated Experiments

Remote sensing datasets have already been used in the current literature. Synthetic data were first considered in [36] and involved geographic images acquired by remote sensing instruments, where the synthesized remote sensing image has a size of 145 × 145 × 724 pixels. The next dataset employed is the Hydice urban dataset [37], with a raw image of 307 × 307 × 210 pixels. In this dataset, the image presenting the entire city is affected by water absorption, the atmosphere, deadlines, and stripes.

To comprehensively evaluate the proposed algorithm, it is challenged against four different image reconstruction methods, i.e., using TV and low-rank matrix factorization tensor completion (LRTV) [35], adaptive weighted tensor completion (AWTC) [16], decomposable nonlocal tensor dictionary learning (TDL) [15] and stochastic gradient descent TT completion (TTSGD) [38, 39]. In the subsequent experiments, the parameters of all competitor techniques are manually adjusted based on the default strategy. In addition, the spectral bands of all remote sensing image sets are normalized to [0, 1] and stretched to the raw level after restoration to facilitate visualization and numerical calculation.

4.1.1. Experiments on the Synthetic Datasets

To simulate noisy remote sensing images and quantitatively highlight the image recovery capabilities of each method, four noise cases are artificially created in which the raw synthetic remote sensing data are polluted with several noise types. In any of these cases, the noise intensity of disparate frequency bands is the same. The added noise on each frequency band has the same distributed zero-mean Gaussian noise, while variance and data loss rate vary. Specifically, the values of variance and data loss rate are, for Case 1, 0.06 and 0.2, for Case 2, 0.03 and 0.2, for Case 3, 0.06 and 0.4, and for Case 4, 0.03 and 0.4, respectively.

The average structural similarity index (MSSIM) and the average peak signal-to-noise ratio (MPSNR) are employed to objectively evaluate the repair quality for understanding and analyzing the overall property of the disparate methods:where and represent the PSNR and SSIM values of the i-th band, respectively.

Figure 2 highlights that all methods have eliminated noise to a certain extent. The proposed method performs best, as it retains the original image’s basic structure while effectively removing noise nuisances. Specifically, ALS obtains the best visual effect, while LRTV, AWTC, TTSGD, and TDL successfully reconstruct the missing stripes. However, there are still noticeable streaks in the results, as some details are smoothed out from the restored image utilizing TTSGD and TDL. Although LRTV produces some excessively smooth details, the methods based on regularization performs best in preserving details and reconstructing fringes. Combining the adaptive nuclear norm regularization algorithm and the Bayesian probability tensor modeling technology affords a decisive function in remote sensing image reconstruction.

Table 1 records three indicators of the restoration results of all competitor methods. The two methods proposed have achieved excellent results in all cases, with ALM being better than ALS, mainly because the adaptive nuclear norm contributes more to the ALM architecture. It is also evident that the improvement afforded by the nuclear norm regularization in the ALS method is not adequate. Additionally, when the missed detection rate is reduced, the AWTC method is more effective than TDL and TTSGD. The results of ALS and ALM are better than LRTL, mainly because the performance of the adaptive nuclear norm normalization is better than the traditional TV regularization. Considering the quantitative evaluation utilizing different indicators on the simulated restoration results, AWTC has the best MPSNR and SLM evaluation, but ALM manages the best MSSIM value.

4.1.2. Experiments on Hydice Urban Datasets

In these experiments, we select the tensor of the decomposed by Bayesian probability RANK [25, 10, 25], [25, 10, 15], [15, 10, 25], and [15, 10, 15], and the adaptive nuclear norm regularize parameter is defined as and for data missing percentages of 25% and 50%. Also, several methods are employed, such as LRTL, AWTC, TDL, and TTSGD.

Table 2 lists the quantitative evaluation results of a random loss of remote sensing images under different completion methods. The table highlights that ALM and ALS have considerable advantages over other comparison methods. The experimental results reveal that under the condition of missing information, the adaptive nuclear norm regularization can effectively enhance the property of low-rank tensor models. The AWTC method achieves better results with a lower miss rate, while the TDL method operates oppositely.

Figure 3 illustrates the composite figure after applying remote sensing image reconstruction when 50% of the data is missing. The figure shows that TDL and TTSGD do not recover some of the missing data, and LRTL and AWTC produce ambiguous results. However, ALM and ALS produce better results. In summary, ALM and ALS perform best for remote sensing reconstruction, demonstrating the advantages of the Bayesian probability tensor decomposition combined with an adaptive nuclear norm adjustment.

4.2. Real-Data Experiments

This section involves experiments utilizing remote sensing camera data from a lake area under actual conditions, aiming to verify the reconstruction capability of the proposed method in real-world imagery. The simulated remote sensing image acquisition device of Figure 4 is used to obtain the corresponding image data, which are further processed to simulate the remote sensing image obtained by the remote sensing device in reality [40]. The gray value of each image band is normalized to [0, 1] before restoration. The tested image is of size of 300 × 300 × 6 × 6, that is, the space size is 300 × 300, the number of spectra is 6, and the time node is 6. We perform cloud pollution on the experimental images to simulate the effects of mixed noise during the remote sensing imagery acquisition.

This trial considers only the recovered images of the models for visual comparison as all competitor methods have similar performance in these actual data experiments. Figure 5 presents the remote sensing images of the time nodes after reconstruction utilizing various methods. The method proposed in this article manages better results compared to the competitor techniques, with ALS achieving the best visual effects. Despite LRTV, TTSGD, AWTC, and TDL can reconstruct the lacking stripes, some noticeable streaks remain. Considering TTSGD and TDL, some details are smoothed out from the enlarged image, and although LRTV produces some overly smooth details, regularization-based methods perform best in maintaining detail and reconstructing fringes.

4.3. Discussion

This subsection considers the computational burden, initialization, convergence, and parameter selection of ALM. Given that ALS and ALM have the same performance, analyzing the latter is neglected to optimize readability and maintain the paper at a reasonable length.

Convergence: Figure 6 presents the interplay between MPSNR as the iteration number in the disparate emulated experiments increases. As the number of iterations increases, the ALM converges for an increasing MPSNR, which at some point stabilizes. Hence, in all trials, 100 iterations are considered.

Sensitivity analysis of Parameter λ: Figure 7 highlights the relationship between MPSNR and parameter λ for real-data experiments. From the latter figure, when parameter λ starts to increase, the value of MPSNR increases and reaches a peak value at 5e − 6, indicating the efficiency of nuclear norm regularization. Nevertheless, the MPSNR value drops due to excessive regularization smoothness after the peak point.

The validity of the Rank Constraint: the Bayesian probability tensor factorization is employed to encode the low-rank prior during the ALM optimization process. Figure 8 shows the MSSIM and MPSNR values of the adaptive nuclear norm solver with disparate rank-constrained values. The results indicate that the values of MSSIM and MPSNR have a specific relationship with the estimated rank value, i.e., both values initially increase and afterward decrease because the estimated rank value increases.

Time Cost: the time cost of disparate methods corresponding to disparate simulation methods is reported in Table 3, demonstrating that the proposed algorithm (ALM) has a lower time complexity than competitor methods.

5. Conclusion

This article proposes a new tensor completion method, where the Bayesian probability tensor decomposition method is used for remote sensing image reconstruction to detect temporal, spectral, and spatial information. In our scheme, we adopt the adaptive kernel-mode regularization and Bayesian probability tensor decomposition for the spatial smoothness of remote sensing images. We demonstrate that the proposed method is more appealing than commonly used methods in visual comparison and quantitative evaluation through a series of simulations and real-data experiments. The conclusions derived from our research are as follows:(i)The suggested method is significantly more robust than the tensor completion method based on traditional tensor transformation. This is proven through our experiments utilizing real or synthetic data. The adaptive nuclear norm uses the gradient similarity between all remote sensing image bands to characterize the smooth structure further, thereby improving the performance and efficiency of the suggested method.(ii)The proposed combination of adaptive nuclear norm regularization and tensor decomposition manages appealing stability and can solve the inherent problems of traditional tensor reconstruction methods, such as high-rank selection sensitivity and low computational efficiency.

Future research shall enhance Bayesian probability tensor decomposition efficiency and extend our method to utilize the gradient similarity between all remote sensing image bands to characterize the smooth structure further.

Abbreviations

:Higher-order tensor
:Latent tensor
:Gaussian distribution
:Accuracy
, :Hyperparameters
:Degree of freedom of the Wishart distribution
:Vertical image
:Mixed noise items
:Gaussian noise term
:Sparse noise term
:Approximate tensor
:Potential factor
:Low-rank tensor
:Weighted nuclear norm
:Weighted N-dimensional difference operator
:Penalty parameter
:Lagrange multipliers
:Elementwise square
:Threshold shrinkage processing
:Elementwise multiplication.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Authors’ Contributions

LJL investigated the study and reviewed and edited the article. YY and SS conceptualized and visualized the study, and reviewed and edited the manuscript. PL and XHQ visualized and investigated the study. CL supervised the study. SL proposed the methodology, was responsible for data curation, performed formal analysis, and wrote the original draft.

Acknowledgments

This research was financially supported by Changchun Science and Technology Bureau Key Research and Development Projects (Grant no. 21ZGG08).