Thick Cloud Removal in Multi-Temporal Remote Sensing Images via Frequency Spectrum-Modulated Tensor Completion

: Clouds often contaminate remote sensing images, which leads to missing land feature information and subsequent application degradation. Low-rank tensor completion has shown great potential in the reconstruction of multi-temporal remote sensing images. However, existing methods ignore different low-rank properties in the spatial and temporal dimensions, such that they cannot utilize spatial and temporal information adequately. In this paper, we propose a new frequency spectrum-modulated tensor completion method (FMTC). First, remote sensing images are rearranged as third-order spatial–temporal tensors for each band. Then, Fourier transform (FT) is introduced in the temporal dimension of the rearranged tensor to generate a spatial–frequential tensor. In view of the fact that land features represent low-frequency components and ﬁckle clouds represent high-frequency components in the time domain, we chose adaptive weights for the completion of different low-rank spatial matrixes, according to the frequency spectrum. Then, Invert Fourier Transform (IFT) was implemented. Through this method, the joint low-rank spatial–temporal constraint was achieved. The simulated data experiments demonstrate that FMTC is applicable on different land-cover types and different missing sizes. With real data experiments, we have validated the effectiveness and stability of FMTC for time-series remote sensing image reconstruction. Compared with other algorithms, the performance of FMTC is better in quantitative and qualitative terms, especially when considering the spectral accuracy and temporal continuity.


Introduction
Satellite remote sensing images have been widely used in many fields, such as geography, ecology and environment monitoring [1,2], and have become one of the most important means of obtaining information on the Earth's surface.However, satellite remote sensing images are prone to contamination by clouds, which causes great difficulties in target detection [3,4], identification, feature classification [5,6] and other applications [7,8].For this reason, remote sensing image inpainting has become an active research area [9].Many methods have been proposed to deal with the reconstruction of missing areas due to cloud contamination.Depending on the information used, these methods can be classified into three categories: spatial-based, spectral-based and temporal-based.
The principle of spatial-based methods is to fill in the missing areas in remote sensing images by propagating surrounding similar structures and texture information.These methods include image interpolation [10,11], propagated methods [12], compressive sensing methods [13], group-based methods [14] and variation-based methods [15,16].However, the spatial-based methods can only be used to reconstruct small missing areas.When the area contaminated by clouds is large, these methods are not applicable.
The spectral-based methods aim to reconstruct the missing areas of remote sensing images using the correlation between different bands or different sensors' images, especially for multispectral images.Among these methods, Aqua MODIS band 6 inpainting is the most typical.For instance, Rakwatin et al. [17] reconstructed the missing data in MODIS band 6 based on the correlation between band 6 and band 7. Roy et al. [18] used information observed by MODIS to predict Landsat ETM images.Since the infrared band can be used to obtain information on the land under thin clouds, Li et al. [19] used the data of the infrared band to reconstruct the cloud-contaminated area of a visual image.In addition, homomorphic filtering [20] and haze-optimized transform (HOT) [21] have been used to deal with thin cloud-contaminated image reconstruction.Since almost no optical band can penetrate the thick clouds, these methods work less well when images are contaminated by thick clouds.
The temporal-based methods use information from time-series remote sensing images in the same regions to reconstruct missing areas [22].Melgani et al. [23] reconstructed the missing regions via an unsupervised contextual prediction process, which uses local spectral-temporal relations at different times.Similarly, Zhang et al. [24] and Lin et al. [25] reconstructed contaminated areas via the correlation of information from other temporal images.In addition, the multi-temporal dictionary learning algorithm was also used for image reconstruction [26].These methods only consider information from the temporal dimension, and fail to take advantage of information in other dimensions, such that the reconstructed results will be limited by the availability of cloud-free areas.
To better utilize the information from the spatial, spectral and temporal dimensions for image reconstruction, many studies have used the low-rank tensor completion method [27].For example, considering the contribution of the spatial, spectral and temporal information in each dimension, Ng et al. [28] proposed an adaptive weighted tensor completion method.Ji et al. [29] proposed a non-local low-rank tensor completion method, which introduced the non-convex approximation of tensor rank and rearranged the fourth-order tensor in groups.Chen et al. [30] considered spatial-spectral total variance regularized low-rank sparsity decomposition for cloud and shadow removal.Chu et al. [31] designed a novel spatialtemporal adaptive tensor completion method to reconstruct the data of cloud-prone regions.Duan et al. [32] proposed a tensor optimization model based on temporal smoothness and sparse representation (TSSTO) for thick cloud and cloud shadow removal in remote sensing images.Lin et al. [33] proposed a robust thick cloud/shadow removal (RTCR) method using coupled tensor factorization to meet the problem arising from an inaccurate mask, and thus reconstruct the multi-temporal information.Low-rank tensor approximation (LRTA) is an emerging technique [34], having gained much attention in the hyperspectral imagery (HSI) restoration community.Liu et al. [35] proposed a multigraph-based lowrank tensor-approximation method for HSI restoration, which integrated geometry-related information into LRTA to constrain a smooth solution of the restored HSIs.
Almost all existing tensor completion methods require the following two steps.First, singular-value threshold decomposition is performed for each mode of the tensor.Second, all tensors are weight-summed.However, for multi-temporal remote sensing images, the low-rank properties in the spatial dimension are different from those in the temporal dimension, which makes the selection of weight difficult.Therefore, in existing tensor completion methods used for remote sensing images, the temporal and spatial information is not utilized reasonably and effectively, which results in reconstructed images that are unclear or have spectral inaccuracy.In order to effectively utilize different low-rank properties in spatial and temporal dimensions, here, we propose a new frequency spectrummodulated low-rank tensor completion method (FMTC).
FMTC treats time-series remote sensing images as the fourth-order tensor, and considers different low-rank processing approaches in the spatial and temporal dimensions.For the frequency spectrum of time-series remote sensing images, the temporal continuity, correlation and periodicity of land features are reflected in the low-frequency part, while clouds and noise influence the high-frequency part.On this basis, remote sensing images are rearranged as third-order spatial-temporal tensors in each band.The Fourier transform is introduced in the temporal dimension of the spatial-temporal tensor, which is converted to the spatial-frequential tensor.Low-pass filtering and the adaptive weights for each spatial matrix are performed via low-rank processing.After that, IFT is implemented.Through this method, joint spatial-temporal low-rank information is derived for reconstructing the missing areas.The main contributions of this paper can be listed as follows:

•
Different orthogonal decompositions are performed on spatial and temporal dimensions of the tensor.The spatial-temporal tensor for each band is transformed to a spatial-frequential tensor by FT.Singular-value decomposition is performed in lowrank matrix completion in the spatial dimension.Through the frequency spectrummodulating spatial matrix, joint spatial-temporal low-rank information is achieved, and the effects of different spatial-temporal low-rank properties are avoided.Meanwhile, using the property of conjugated symmetry of FT can also reduce the computation cost during the iteration.

•
Gaussian low-pass filtering is applied in the frequency spectrum, and spatial low-rank adaptive weights are calculated according to the frequency characteristics of the time domain.Thus, the difficulty in selecting appropriate weights is solved.This scheme can maintain the low-frequency land features and weaken the high-frequency noise caused by clouds.
The rest of this paper is organized as follows: Section 2 introduces the basic idea, the model of tensor completion and the algorithm of this paper.Section 3 shows the results and analysis of simulated and real data experiments.Section 4 gives the conclusion.

Spatial-Temporal Low-Rank Tensor Rearrangement
For multi-temporal remote sensing images obtained by satellites with high revisiting frequency, land features slowly change or are invariant over a period of time.Clouds break the temporal continuity of land features.Meanwhile, information on a region contaminated by thick clouds is lost in all bands.Therefore, more effective information for the reconstruction of images is given in the temporal dimension.
The original multi-temporal remote sensing images can be represented as Y ∈ R m×n×b×t , where m and n are the spatial sizes of the image, b is the number of bands, and t is the number of time-series images.In order to make use of the information in the temporal dimension, Y is rearranged, as shown in Figure 1.The data of each band are expressed as a spatial-temporal third-order tensor, and the number of tensors is b.In this paper, we handle the spatial-temporal tensor for each band separately.For convenience, we denote the spatial-temporal tensor as X ∈ R m×n×t .The initial spatial-temporal tensor is denoted as X 0 .Ω ∈ R m×n×t is the cloud mask tensor of X 0 , where Ω ijk = 1 marks clear pixels and Ω ijk = 0 marks the pixels contaminated by clouds, and (i, j, k) are the discrete indices of (m, n, t).
Due to the correlation and continuity of remote sensing data in the spatial and temporal dimensions, the rearranged spatial-temporal tensor X ∈ R m×n×t has low-rank properties.The following low-rank tensor completion model [36] can be applied to reconstruct images, min where X, T ∈ R m×n×t , T Ω is the spatial-temporal tensor, including pixels that are not obscured by clouds.Instead of rank(X) of the trace norm, the optimization problem is changed to the following form [27]: where X * is the trace norm of X.The unconstrained version of the problem [27] can be written as min where ρ is a pre-set constant, and F is the Frobenius norm of the tensor.In the form of the solution of Equation (3), the "shrinkage" operator D ω (X) [37] is introduced into the low-rank processing of X, where ω is the weight used in the calculation of the trace norm for each mode.It is difficult for the spatial-temporal tensor to select appropriate weights for different spatial-temporal low-rank properties.FMTC is proposed to solve the problem.Due to the correlation and continuity of remote sensing data in the spatial and temporal dimensions, the rearranged spatial-temporal tensor ∈ ℝ × × has low-rank properties.The following low-rank tensor completion model [36] can be applied to reconstruct images, min : rank( ) . .: = , where , ∈ ℝ × × , is the spatial-temporal tensor, including pixels that are not obscured by clouds.Instead of rank( ) of the trace norm, the optimization problem is changed to the following form [27]: where ‖ ‖ * is the trace norm of .The unconstrained version of the problem [27] can be written as where is a pre-set constant, and ‖ − ‖ is the Frobenius norm of the tensor.In the form of the solution of Equation (3), the "shrinkage" operator ( ) [37] is introduced into the low-rank processing of , where is the weight used in the calculation of the trace norm for each mode.It is difficult for the spatial-temporal tensor to select appropriate weights for different spatial-temporal low-rank properties.FMTC is proposed to solve the problem.

Frequency Spectrum-Modulated Tensor Completion
Fourier transform is a method for projecting time-domain signals onto a set of orthogonal trigonometric bases.It is especially suitable for periodic data decomposition and processing.FMTC uses the frequency spectrum in the temporal dimension to adaptively determine the low-rank weight of each spatial slice matrix, which preserves the low-frequency components of the land features in the temporal dimension and limits the highfrequency noise caused by clouds.Then, we can achieve the joint spatial-temporal low-

Frequency Spectrum-Modulated Tensor Completion
Fourier transform is a method for projecting time-domain signals onto a set of orthogonal trigonometric bases.It is especially suitable for periodic data decomposition and processing.FMTC uses the frequency spectrum in the temporal dimension to adaptively determine the low-rank weight of each spatial slice matrix, which preserves the low-frequency components of the land features in the temporal dimension and limits the high-frequency noise caused by clouds.Then, we can achieve the joint spatial-temporal low-rank constraint.
First, we perform the Fourier transform on the temporal dimension of the rearranged tensor X ∈ R m×n×t to generate a spatial-frequential tensor, where [] denotes the transform length, which is the default value in matlab.
Then, for the spatial slice matrix Xi ∈ R m×n of the spatial-frequential tensor X, i = 1, • • • , t, the singular-value decomposition is performed: where U, S, V are the matrices of Xi 's singular-value decomposition [38].
Examples of a rearranged spatial-temporal tensor and the spatial-frequential tensor are shown in Figure 2. The low-frequency part represents slowly changing or constant land features, shown in the first two images in Figure 2b, and the high-frequency part represents significantly changed land features, clouds and noise, as shown in the last two images in Figure 2b. Figure 3b shows the time-series scatter plot of the near-infrared band from 2003 to 2008 for one pixel of Figure 3a, where the value 0 represents contamination by clouds or data invalidation.After removing the influence of clouds, the periodicity of land features that changes with seasons can be seen in the long-term image sequences.Figure 4 shows the Fourier transform spectrum curve of Figure 3.The maximum value represents zero frequency.The two large values denoted by red dots in the low-frequency part are caused by the periodic variation of land features.Noise and clouds appear in the high-frequency part.
where , , are the matrices of 's singular-value decomposition [38].Examples of a rearranged spatial-temporal tensor and the spatial-frequential tensor are shown in Figure 2. The low-frequency part represents slowly changing or constant land features, shown in the first two images in Figure 2b, and the high-frequency part represents significantly changed land features, clouds and noise, as shown in the last two images in Figure 2b. Figure 3b shows the time-series scatter plot of the near-infrared band from 2003 to 2008 for one pixel of Figure 3a, where the value 0 represents contamination by clouds or data invalidation.After removing the influence of clouds, the periodicity of land features that changes with seasons can be seen in the long-term image sequences.where , , are the matrices of 's singular-value decomposition [38].Examples of a rearranged spatial-temporal tensor and the spatial-frequential tensor are shown in Figure 2. The low-frequency part represents slowly changing or constant land features, shown in the first two images in Figure 2b, and the high-frequency part represents significantly changed land features, clouds and noise, as shown in the last two images in Figure 2b. Figure 3b shows the time-series scatter plot of the near-infrared band from 2003 to 2008 for one pixel of Figure 3a, where the value 0 represents contamination by clouds or data invalidation.After removing the influence of clouds, the periodicity of land features that changes with seasons can be seen in the long-term image sequences.According to frequency spectrum characteristics of the spatial-frequential tensor, we proposed the following operations for .First, in order to maintain the low-frequency land features and weaken the high-frequency noise caused by clouds, a Gaussian lowpass filter is applied to in the frequency spectrum.Second, spatial adaptive weights According to frequency spectrum characteristics of the spatial-frequential tensor, we proposed the following operations for X.First, in order to maintain the low-frequency land features and weaken the high-frequency noise caused by clouds, a Gaussian lowpass filter is applied to X in the frequency spectrum.Second, spatial adaptive weights determined by the frequency spectrum are performed to achieve the joint spatial-temporal low-rank constraint.
According to the properties of Fourier transform, there is a conjugate symmetry in the frequency domain of X.The Gaussian filter is designed as follows: , and σ is a pre-set constant.We define each spatial matrix of X as Xi , i = 1, • • • , t. Due to Xi representing different frequency parts, when the same threshold is applied to perform the spatial low-rank processing of Xi , it will cause the loss of land feature information.In order to maintain land feature information in the image reconstruction, we choose the adaptive weight ω i for Xi according to the importance of the information from X1 to Xt .
Xi is defined as the normalized mean value k i as follows: where m i denotes the mean value of Xi .A larger k i means more important information is contained in Xi , and a smaller ω i should be used in the low-rank processing.The adaptive weight ω i of each spatial matrix is defined as follows: where ρ is a pre-set constant.

Model Optimization
We define the augmented Lagrangian equation as follows: where M, B and X are updated by the alternating direction method of multipliers (ADMM) [39].M has a close-form solution: where k is the iteration time.ifft(•) is the Inverse Fast Fourier Transform.The adaptive weight ω i is determined by Equation (8).For X, it is updated as follows: For B, it is updated as follows: In this paper, the initial ρ is set to 10 −4 .Later, to accelerate convergence, ρ will be iteratively increased by ρ k+1 = tρ k , t ∈ [1.1, 1.2].The solution process of FMTC is shown in Algorithm 1.

Experimental Data
In this part, Landsat Collection 1 L1TP data from 2003 to 2018 in Beijing have been selected as the experiment data, as shown in Table 1.The experiment data include five simulated datasets and two real datasets based on different land-cover types and cloud areas, respectively.In order to demonstrate the effectiveness of FMTC, five algorithms are selected for comparison, including HaLRTC [27], AWTC [28], NL-LRTC [29], TVLRSD [30] and ST-Tensor [31].

Evaluation Indicators
In the simulated experiments, the Peak Signal-to-Noise Ratio (PSNR) [40], Structural Similarity (SSIM) [41] and Spectral Angle Mapper (SAM) [42] are chosen to evaluate the reconstructed images from the spatial and spectral perspectives.
The PSNR is calculated as follows: where MAX is the maximum pixel value of the original image, and I(i, j) and K(i, j) denote the (i, j)th pixel values of the original image and the reconstructed image, respectively.m and n are the spatial sizes of those images.The SSIM is calculated as follows: where µ x and µ y represent the mean values of the original image and the reconstructed image, respectively.σ x and σ y represent the standard deviations.σ xy represents the covariance.
The SAM is calculated as follows: where (x 1 , • • • , x n ) and (y 1 , • • • , y n ) denote the original image and reconstructed image of the spectral vectors.
In the real data experiments, the above evaluation indicators cannot be used because there was no reference image.Therefore, the information entropy (IE) [43] and the average gradient (AG) [44] are used to evaluate the reconstructed image.In addition, the effectiveness of the algorithm can be evaluated by the consistency of spatial images, the clarity of feature targets and the continuity of the time-series curve.
Information entropy of the reconstructed image is calculated as follows: where p t is the probability of value t.G is the maximum pixel value of the reconstructed image.
The average gradient of the reconstructed image is calculated as follows:

Simulated Data Experiments
In the simulated experiments, the effect of clouds is simulated by masking a random area on the cloud-free images.In order to demonstrate the effectiveness of FMTC on different land-cover types and missing areas, two simulated experiments were conducted.One is based on datasets 1-4 with four land-cover types, the other is based on dataset 5 with different sizes of missing areas.The information of the datasets is shown in Table 2.  5a, 7a, 9a and 11a are true-color composition images (red-band 3, green-band 2, and blue-band 1).The simulated data have been generated from random missing areas where the pixel values were set to 0, as shown in Figures 5b, 7b According to the enhanced detail shown in Figures 5, 7, 9 and 11, we can see that the reconstructed results of HaLRTC and AWTC are blurry for four land-cover types because these algorithms only utilize the low-rank information derived from experimental data.The results of NL-LRTC are better than those of the first two algorithms, but still a bit fuzzy, which is due to the inappropriate choice of weights for the different modes of the tensor.TVLRSD can reconstruct the basic structure and details of the images, but the reconstructed details are still not clear.This may be because the matrix decomposition denoising the image and removing clouds makes images smooth and vague, but the details are not reconstructed.Although ST-Tensor maintains much of the detail of the images, the color of the reconstructed areas is inconsistent with that of the surroundings, due to the inappropriate choice of weights in the rearranged spatial-temporal tensor.In contrast, FMTC is able to ensure spatial detail as well as spectral consistency.
fuzzy, which is due to the inappropriate choice of weights for the different modes of the tensor.TVLRSD can reconstruct the basic structure and details of the images, but the reconstructed details are still not clear.This may be because the matrix decomposition denoising the image and removing clouds makes images smooth and vague, but the details are not reconstructed.Although ST-Tensor maintains much of the detail of the images, the color of the reconstructed areas is inconsistent with that of the surroundings, due to the inappropriate choice of weights in the rearranged spatial-temporal tensor.In contrast, FMTC is able to ensure spatial detail as well as spectral consistency.A quantitative comparison for the six algorithms applied to four land-cover types is shown in Table 3. HaLRTC, AWTC and NL-LRTC obtained poor results for all quantitative indicators.The results of TVLRSD and ST-Tensor are better than those of the above three algorithms.Generally speaking, FMTC can achieve much better results than the other five algorithms.Although the PSNR and SSIM of ST-Tensor seem to perform a little better than FMTC on impervious land, the SAM of ST-Tensor is significantly effective than FMTC.In addition, FMTC shows outstanding performance on the three other land-cover types for each indicator.It can be seen that the results of the quantitative evaluation are consistent with visual examination.It is worth mentioning that FMTC is far less time-consuming than the other five algorithms.A quantitative comparison for the six algorithms applied to four land-cover types is shown in Table 3. HaLRTC, AWTC and NL-LRTC obtained poor results for all quantitative indicators.The results of TVLRSD and ST-Tensor are better than those of the above three algorithms.Generally speaking, FMTC can achieve much better results than the other five algorithms.Although the PSNR and SSIM of ST-Tensor seem to perform a little better than FMTC on impervious land, the SAM of ST-Tensor is significantly less effective than FMTC.In addition, FMTC shows outstanding performance on the three other landcover types for each indicator.It can be seen that the results of the quantitative evaluation  A quantitative comparison for the six algorithms applied to four land-cover types is shown in Table 3. HaLRTC, AWTC and NL-LRTC obtained poor results for all quantitative indicators.The results of TVLRSD and ST-Tensor are better than those of the above three algorithms.Generally speaking, FMTC can achieve much better results than the other five algorithms.Although the PSNR and SSIM of ST-Tensor seem to perform a little better than FMTC on impervious land, the SAM of ST-Tensor is significantly less effective than FMTC.In addition, FMTC shows outstanding performance on the three other landcover types for each indicator.It can be seen that the results of the quantitative evaluation are consistent with visual examination.It is worth mentioning that FMTC is far less time-     In this simulated experiment, dataset 5 was used to demonstrate the effectiveness of FMTC on different missing sizes.The proportions of the missing area are 6.01%, 19.26% and 32.48%, respectively.Figures 13-15 are visual presentations of the results for three missing sizes, respectively.Figures 13a, 14a and 15a show true-color composition images (red-band 4, green-band 3, and blue-band 2).Simulated cloudy images are shown in Figures 13b, 14b and 15b.Figures 13c, 14c and 15c show the reconstructed images.Figures 13d-f, 14d-f and 15d-f show, in enhanced detail, the areas in the red boxes from Figures 13a-c, 14a-c and 15a-c, respectively.
Visually, the reconstructed results are consistent with the surroundings, and the proportion of the missing size reaches 32.48%.According to the areas with enhanced detail, the results of reconstruction approach the original images.In Figure 13d, a white road is depicted.In Figure 13f, the white road is reconstructed completely, owing to the timeseries information utilized by FMTC.The above results demonstrate the effectiveness of using FMTC with different missing area sizes.A quantitative comparison of the six algorithms with different missing sizes is shown in Table 4.As the proportion of missing size increases, the PSNR and SSIM of the first three algorithms decrease significantly, while the SAM increases.In contrast, the latter   A quantitative comparison of the six algorithms with different missing sizes is shown in Table 4.As the proportion of missing size increases, the PSNR and SSIM of the first three algorithms decrease significantly, while the SAM increases.In contrast, the latter three algorithms are less sensitive to changes in missing sizes.Among all the algorithms, Visually, the reconstructed results are consistent with the surroundings, and the proportion of the missing size reaches 32.48%.According to the areas with enhanced detail, the results of reconstruction approach the original images.In Figure 13d, a white road is depicted.In Figure 13f, the white road is reconstructed completely, owing to the time-series information utilized by FMTC.The above results demonstrate the effectiveness of using FMTC with different missing area sizes.
A quantitative comparison of the six algorithms with different missing sizes is shown in Table 4.As the proportion of missing size increases, the PSNR and SSIM of the first three algorithms decrease significantly, while the SAM increases.In contrast, the latter three algorithms are less sensitive to changes in missing sizes.Among all the algorithms, FMTC performs the best in most of the evaluation indicators.For different missing sizes, FMTC can obtain satisfactory results.

Real Data Experiments
In experiments with real data, two datasets were selected, as shown in Table 5.For the real data, we used a Landsat quality assessment to mask clouds and cloud shadows.In order to demonstrate the effectiveness of applying FMTC on different data sources, two real data experiments were conducted.In this part, dataset 6 from Landsat-5 for the first real data experiment was selected.The main land cover type is impervious, and the others are vegetation and soil.Figure 16a shows the true-color composition image.Figure 16b shows the masked image.Visual comparisons of FMTC with other algorithms are shown in Figure 16c-h.Figure 17a-h show enhanced details of the areas in the red boxes.
The results of this experiment are similar to those of the simulated data experiment.The reconstructed results of HaLRTC and AWTC are not satisfactory and cannot meet the basic demands of reconstruction.The reconstructed details of TVLRSD and NL-LRTC are not clear, and the result of ST-Tensor is inconsistent with its surroundings.FMTC can not only clarify the reconstructed details, but it also ensures the consistency of the spectrum.
A quantitative comparison is shown in Table 6.The IE and AG of FMTC are better than those of other algorithms, indicating that the reconstructed details are richer and clearer.Figure 18 shows the time-series curves of six algorithms applied to one pixel from 2004 to 2006.The value 0 represents that the pixel was contaminated by cloud on that day.The curves of HaLRTC, AWTC and NL-LRTC display values of 0, leading to unsatisfactory results.The curves of TVLRSD and ST-Tensor are better than those of the first three algorithms, while they still show abrupt changes for the cloudy pixels.The curve of FMTC looks more continuous than those of all other algorithms.It can thus better characterize the temporal continuity and periodicity of land features.A quantitative comparison is shown in Table 6.The IE and AG of FMTC are better than those of other algorithms, indicating that the reconstructed details are richer and clearer.Figure 18 shows the time-series curves of six algorithms applied to one pixel from  A quantitative comparison is shown in Table 6.The IE and AG of FMTC are better than those of other algorithms, indicating that the reconstructed details are richer and clearer.Figure 18 shows the time-series curves of six algorithms applied to one pixel from   In this part, dataset 7 from Landsat-8 was selected for use in the second real data experiment.The main land cover type is vegetation, and the others are impervious and soil.Figure 19a is the original image.Figure 19b is the masked image.Visual comparisons of FMTC with other algorithms are shown in Figure 19c-h.Figure 20a-h show enhanced details of the areas in the red boxes.
In terms of the visual effect, the results of this experiment are similar to those from real data experiment 1.As regards the areas of enhanced detail, TVLRSD is not clear enough.The white road, shown in the red box in Figure 20g,h, is not obvious in the reconstructed result of ST-Tensor, while the white road is reconstructed completely by FMTC.A quantitative evaluation is shown in Table 7.The time-series curves of the six algorithms applied to one pixel from 2015 to 2017 are shown in Figure 21.FMTC achieves the best performance, with excellent reconstructions from different data sources.

Real Data Experiment 2
In this part, dataset 7 from Landsat-8 was selected for use in the second real data experiment.The main land cover type is vegetation, and the others are impervious and soil.Figure 19a is the original image.Figure 19b is the masked image.Visual comparisons of FMTC with other algorithms are shown in Figure 19c-h.Figure 20a-h

Real Data Experiment 2
In this part, dataset 7 from Landsat-8 was selected for use in the second real data experiment.The main land cover type is vegetation, and the others are impervious and soil.Figure 19a is the original image.Figure 19b is the masked image.Visual comparisons of FMTC with other algorithms are shown in Figure 19c-h.Figure 20a-h show enhanced details of the areas in the red boxes.
In terms of the visual effect, the results of this experiment are similar to those from real data experiment 1.As regards the areas of enhanced detail, TVLRSD is not clear enough.The white road, shown in the red box in Figure 20g,h, is not obvious in the reconstructed result of ST-Tensor, while the white road is reconstructed completely by FMTC.A quantitative evaluation is shown in Table 7.The time-series curves of the six algorithms applied to one pixel from 2015 to 2017 are shown in Figure 21.FMTC achieves the best performance, with excellent reconstructions from different data sources.

Conclusions
In this paper, we propose a frequency spectrum-modulated tensor completion method for multi-temporal remote sensing image inpainting.We rearrange the original four-dimensional tensor to establish a series of spatial-temporal tensors with low-rank properties.Fourier transform is introduced in the temporal dimension to generate the spatial-frequential tensor.The spatial adaptive weights are determined by the spectral properties of the time-domain signal, which is filtered using the Gaussian low-pass filtering to In terms of the visual effect, the results of this experiment are similar to those from real data experiment 1.As regards the areas of enhanced detail, TVLRSD is not clear enough.The white road, shown in the red box in Figure 20g,h, is not obvious in the reconstructed result of ST-Tensor, while the white road is reconstructed completely by FMTC.A quantitative evaluation is shown in Table 7.The time-series curves of the six algorithms applied to one pixel from 2015 to 2017 are shown in Figure 21.FMTC achieves the best performance, with excellent reconstructions from different data sources.

Conclusions
In this paper, we propose a frequency spectrum-modulated tensor completion method for multi-temporal remote sensing image inpainting.We rearrange the original four-dimensional tensor to establish a series of spatial-temporal tensors with low-rank

Figure 3 . 21 Figure 4 .
Figure 3. (a) Original image (20030327).(b) Scatter plot of the time series of the near-infrared band from 2003 to 2008 for a pixel of (a) in the red box.Figure 3. (a) Original image (20030327).(b) Scatter plot of the time series of the near-infrared band from 2003 to 2008 for a pixel of (a) in the red box.Remote Sens. 2023, 15, x FOR PEER REVIEW 6 of 21

Figure 4 .
Figure 4.The frequency spectrum curve of the above data in Figure 3.

Algorithm 1
Optimization of the FMTC method Input: Original remote sensing data Y, parameter ρ.Initialize: M 0 = B 0 = 0, max iteration times K = 200.for i = 0 to b do Obtain the rearranged spatial-temporal tensor X 0 for each band.for k = 0 to K do Update M by Equation (10); Update X by Equation (11); Update B by Equation (12); end for; repeat until convergence; end for; Output: Reconstructed image data X.
, 9b and 11b.Visual comparisons of FMTC with other algorithms are shown in Figure 5c-h, Figure 7c-h, Figures 9c-h and 11c-h.Figures 6, 8, 10 and 12 show enhanced details of the images.
Figures 13-15 are visual presentations of the results for three missing sizes, respectively.Figures 13a, 14a and 15a show true-color composition images (red-band 4, green-band 3, and blue-band 2).Simulated cloudy images are shown in Figures 13b, 14b and 15b.Figures 13c, 14c and 15c show the reconstructed images.
Figure 13d-f, Figures 14d-f and 15d-f show, in enhanced detail, the areas in the red boxes from Figure 13a-c, Figures 14a-c and 15a-c, respectively.

Figure 13 . 21 Figure 13 .
Figure 13.The simulated experiment with a missing area of 6.01%: (a) original image; (b) simulated cloudy image; (c) reconstructed image; (d-f) enhanced details from red boxes in (a-c).

Figure 14 .
Figure 14.The simulated experiment with missing area of 19.26%:(a) original image; (b) simulated cloudy image; (c) reconstructed image; (d-f) enhanced details from red boxes in (a-c).

Figure 15 .
Figure 15.The simulated experiment with missing area of 32.48%:(a) original image; (b) simulated cloudy image; (c) reconstructed image; (d-f) enhanced details from red boxes in (a-c).

Figure 14 .
Figure 14.The simulated experiment with missing area of 19.26%:(a) original image; (b) simulated cloudy image; (c) reconstructed image; (d-f) enhanced details from red boxes in (a-c).

Figure 15 .
Figure 15.The simulated experiment with missing area of 32.48%:(a) original image; (b) simulated cloudy image; (c) reconstructed image; (d-f) enhanced details from red boxes in (a-c).

Figure 15 .
Figure 15.The simulated experiment with missing area of 32.48%:(a) original image; (b) simulated cloudy image; (c) reconstructed image; (d-f) enhanced details from red boxes in (a-c).

Figure 18 .
Figure 18.The time-series curves of six algorithms applied to one pixel from 2004 to 2006.
show enhanced details of the areas in the red boxes.

Figure 18 .
Figure 18.The time-series curves of six algorithms applied to one pixel from 2004 to 2006.

Figure 21 .
Figure 21.The time-series curves of six algorithms applied to one pixel from 2015 to 2017.

Table 1 .
Experiment data information.

Table 2 .
Information on five datasets.The Experiment Based on Different Land-Cover Types In this part, Landsat-5 datasets for four different regions in Beijing from 2003 to 2011 were selected, as shown in Table 2.The experiment results are shown in Figures 5-12.Figures

Table 3 .
Quantitative comparison of six algorithms applied to four land-cover types.

Table 3 .
Cont.The Experiment with Different Missing SizesIn this simulated experiment, dataset 5 was used to demonstrate the effectiveness of FMTC on different missing sizes.The proportions of the missing area are 6.01%, 19.26% and 32.48%, respectively.

Table 4 .
Quantitative comparison of six algorithms with different missing sizes.

Table 6 .
Quantitative comparison of six algorithms applied to dataset 6.

Table 6 .
Quantitative comparison of six algorithms applied to dataset 6.
Figure 18.The time-series curves of six algorithms applied to one pixel from 2004 to 2006.3.4.2.Real Data Experiment 2

Table 7 .
Quantitative comparisons of six algorithms applied to dataset 7. The time-series curves of six algorithms applied to one pixel from 2015 to 2017.

Table 7 .
Quantitative comparisons of six algorithms applied to dataset 7. Remote Sens. 2023, 15, x FOR PEER REVIEW 18 of 21

Table 7 .
Quantitative comparisons of six algorithms applied to dataset 7.
Figure 21.The time-series curves of six algorithms applied to one pixel from 2015 to 2017.