## 1 Introduction  
Remote sensing datasets provide long records of earth observation, which are essential for earth sciences research. Since the 1960s, remote sensing satellite datasets (e.g., NOAA/AVHRR, MODIS, Landsat) have been widely applied to the study of global environmental change [1; 2; 3; 4]. With the development of remote sensing retrieval algorithms, a series of land surface variables products have been archived, such as the normalized difference vegetation index (NDVI), leaf area index (LAI), gross primary productivity (GPP), and atmospherically corrected surface reflectance [5; 6; 7]. Nonetheless, the time series land surface variables are not always evenly distributed and continuous due to cloud, cloud shadows, snow contamination, low temporal frequency, sun-sensor-surface viewing geometry of different seasons, different sensors and sensor failure, etc. High-precision land cover classification, monitoring, and trend analysis all require continuous fine-scale Remotely Sensed Land Surface Products (RSLSP) [8; 9; 10; 11; 12; 13]. To obtain temporally and spatially continuous RSLSP datasets, spatial-, spectral-, and temporal-based methods and their combinations (e.g., joint spatio-spectral methods and joint spatio-temporal methods) for the time series reconstruction of RSLSP have been proposed [10]. The temporal-based methods have been widely used to remove noises in the time series reconstruction, which can be categorized into five categories [14]: (1) temporal interpolation methods, (2) temporal filters, (3) temporal function-fitting methods, (4) temporal deep learning methods and (5) frequency-based methods. The temporal interpolation methods like Maximum Value Composite (MVC) [15; 16; 17], iterative interpolation for data reconstruction (IDR) [18], data assimilation (DA) [19; 20; 21] and the best index slope extraction (BISE) [22; 23] have been commonly applied to filling time series gaps. However, such algorithms are not suitable for the reconstruction of data with long-term gaps. Due to no consideration of atmosphere interference, the MVC method can lead to the loss of useful temporal information. The IDR method proves prone to overestimating the time series values during the vegetation dormant periods [24]. The DA method has poor performance at the onset of vegetation spring green [20], and the BISE method is bad for long-term decline trends [25]. As for the temporal filter methods, there are the Savitzky-Golay (SG) filter, the Mean-Value iterative (MVI) filter, the Moving Average Filter, and the Changing Weight Filter. The famous Savitzky-Golay (SG) filter method has been widely applied to remove temporal noises from the NDVI, GPP, and reflectance time series [26; 27]. However, the processing effect of the SG filter method is highly dependent on the size of the filtering window size. A large window size can also lead to the existing variations being smoothed out, and a small window size can not effectively reduce the noises in the time series [28; 29]. The MVI filter has negative effects for areas with high interannual changes, and the MA filter will change the width and the height of the curve, for the CW method, the remaining noises in the time series can be problematic. The temporal function fitting method reconstructs the time series by stimulating the time-series variation with parameterized functions. The asymmetric Gaussian (AG) model can extract the phonological parameters during the smoothing process. However, it will lead to a loss of details about vegetation changes if there is a lot of noise in the original time series [30; 31]. The double logistic (DL) technique is effective for revealing the trend of the NDVI time series, but it remains problematic with NDVI in winter [32]. Among frequency-based techniques, the Fourier transform is very effective in processing periodic time-series signals, as it can decompose them into amplitude and phase information, but it has poor performance for irregular time series due to the strict symmetry of the formula [33; 34].