A Continuous Change Tracker Model for Remote Sensing Time Series Reconstruction

: It is hard for current time series reconstruction methods to achieve the balance of high-precision time series reconstruction and explanation of the model mechanism. The goal of this paper is to improve the reconstruction accuracy with a well-explained time series model. Thus, we developed a function-based model, the CCTM (Continuous Change Tracker Model) model, that can achieve high precision in time series reconstruction by tracking the time series variation rate. The goal of this paper is to provide a new solution for high-precision time series reconstruction and related applications. To test the 6-coefﬁcients method. Finally, the CCTM model also has the potential to be applied to change detection, trend analysis, and phenology and seasonal characteristics extractions.


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]. The wavelet transform (WT) method also easily ignores the reasonable high values when there are frequent variations in the time series [35].
The harmonic analysis of time series (HANTS) method was later proposed by Verhoef [36] and has been widely used to reconstruct time-series datasets of the NDVI, LAI, and land surface temperatures (LSTs) [37][38][39]. The method has also been successfully used to predict Landsat values at any given time and has achieved relatively high accuracy in Landsat data reconstruction [9]. However, the HANTS may result in a large deviation from the original time series if there are long-term data gaps [40].
In this paper, a hybrid time series reconstruction method has been proposed. Several time series outlier recognition and gap-filling methods have been applied to fill the longterm gaps. A high-precision model continuous change rate tracking model (CCTM) has been developed for the reconstruction of time series RSLSP. The CCTM model divides the yearly time series into four seasons and reconstructs each seasonal time series with a leastsquares solver. The CCTM model decomposes the land surface variations into four parts by the variation types: uniform variation part, accelerated variation part, decelerated variation part, and short-term periodic variation part. Every part is represented by a designed function. In addition, the SG filter is adapted for LAI, GPP, and surface reflectance time series data to effectively control the degree of smoothing and retain the original characteristics by setting limited iteration times or adding smoothing controlling parameters.
Current time series reconstruction methods do have their advantages, but the longterm gaps, over-fitting, and remaining noise issues are still disturbing [41]. More importantly, none of the methods above can precisely predict the time series values on a daily scale. The objectives of this paper are: (1) to provide solutions to the time series deficiencies of long-term gaps, over-fitting, and remaining noise issues; and (2) to develop a comprehensive continuous time series reconstruction method that can precisely predict time series values.

Data
A series of land surface variable datasets were chosen to evaluate the reconstruction performance of the new model. The NDVI (produced by MOD09Q1), LAI (MCD15A3H), GPP (MOD17A2H), and MSR (MOD09Q1) datasets from 2001 to 2005 in the global region are selected. 183 NDVI images, 183 GPP images, 183 MSR images, and 228 LAI images are collected in the global region. The detailed information is shown in Table 1. Among them, the NDVI and MSR datasets have temporal and spatial resolutions of 8 days and 250 m, the LAI dataset has temporal and spatial resolutions of 4 days and 500 m, and the GPP dataset has temporal and spatial resolutions of 8 days and 500 m.

Time Series Reconstruction Flow
The data processing consists of three parts ( Figure 1). (1) Preprocessing includes setting a mask for the dataset according to the quality control information, outlier recognition, gap-filling, and enhanced SG filtering, which is only necessary for time series datasets that contain noise. If >50% of the points in the yearly time series are invalid, the invalid points are filled with zeros, and gap-filling and enhanced SG filtering are not needed. (2) The reconstruction results are compared in seasonal time series patterns and yearly time-series patterns. In the yearly time series pattern, four models are applied to reconstruct the yearly segment time series, but in the seasonal segment time series, four models are separately used to reconstruct every seasonal segment time series. The MOD09Q1 product provides an estimate of the surface spectral reflectance of bands 1 and 2 at 250 m resolution and corrected for atmospheric conditions such as gasses, aerosols, and Rayleigh scattering.
The NDVI time series dataset is based on the MOD09Q1 product. MCD15A3H, MOD11A2, MOD17A2H, and MOD09Q1 time-series datasets are provided with detailed quality control instructions. The NDVI, LAI, and GPP products were chosen as they are important in reflecting vegetation changes. The MSR product (including near-infrared and red band data) was chosen as it is a representative primary remote sensing product.

Data Preprocessing
Data preprocessing involves two parts: (1) point sampling and quality cont conducted on the Google Earth Engine (GEE) platform; and (2) outlier recognitio filling, and enhanced SG filtering are conducted.  These four-time series reconstruction models include: (i) the CCTM model developed in this paper, (ii) the Harmonic model (an eight-parameter harmonic model), (iii) EXP (an exponential model), and (iv) LN model (a logarithmic model). (3) Accuracy assessment evaluates the reconstruction accuracy between the original time series processed by SG filter and gap-filling methods and the continuous-time series reconstructed by the four models.

Data Preprocessing
Data preprocessing involves two parts: (1) point sampling and quality control are conducted on the Google Earth Engine (GEE) platform; and (2) outlier recognition, gapfilling, and enhanced SG filtering are conducted.

Points Sampling and Quality Control
The six kinds of datasets were processed using the following two steps. First, setting a mask for time series data with the quality control band provides the pixel quality status of each point in the time series. Second, to comprehensively analyze the reconstruction effects on the specific landcover types, the stable land cover type was defined as the land cover which has not been changed from 2001 to 2005, and eight stable land cover types were extracted from the MCD12Q1.006 [42]  Noises in the time series can be reduced by methods, such as the maximum composite method (choosing the highest-quality pixel in a time interval [15]) and the typical time-series filter [36,41,43]. However, there were still spurious points that remained. Since outliers can deviate obviously from a gradually changing time series and always exhibit Noises in the time series can be reduced by methods, such as the maximum composite method (choosing the highest-quality pixel in a time interval [15]) and the typical timeseries filter [36,41,43]. However, there were still spurious points that remained. Since outliers can deviate obviously from a gradually changing time series and always exhibit sudden increases or decreases, a series of algorithms were provided to recognize the spurious points in the time series.
An outlier is defined based on two criteria. The first is that its value is greater/lower than the mean of its moving window plus/minus the cutoff (cutoff is twice the standard deviation of the points in the moving window [43]. In Equation (1), sigma (y(t − c), y(t − c + 1), . . . , y(t), y(t + 1), . . . , y(t + c)) is the standard deviation of the point values in the moving window). It is the serial number of time series points, where y(t) is the value of the tth time series point, c is half of the moving window size, and 2c + 1 is the size of the moving window. y(t) < mean(y(t − c), y(t − c + 1), . . . , y(t), y(t + 1), . . . , y(t + c)) − cuto f f y(t) > mean(y(t − c), y(t − c + 1), . . . , y(t), y(t + 1), . . . , y(t + c)) + cuto f f cuto f = 2 * sigma(y(t − c), y(t − c + 1), . . . , y(t), y(t + 1), . . . , y(t + c)) In Equation (2), N K is the value of the kth point; dN K−1 N K is the inclining or declining slope from the (k − 1)th point N K−1 to the kth point N K ; dN K N K+1 is the inclining or declining slope from the (k − 1)th point N K to the kth point N K+1 ; K is the serial number of the time series point; and ndt is the threshold of the time series points whose interval is n-day, which is defined as 30%/16 × n and equal to 30% when the temporal resolution is 16 days. The second criterion is that both dN K−1 N K and dN K N K+1 are larger than ndt or lower than −ndt.
If one of the above criteria is met, the point in the time series is defined as an outlier. The second criterion is always used in the outlier recognition of NDVI time series and is consistent with the best index slope extraction algorithm (BISE) [23].

Gap Filling
Linear and spline interpolations methods are widely applied to fill the time series gaps [18,26]. For each sampling point, the gaps in the time series were filled using the three gap-fill algorithms presented in this paper. The weighted average method is used to fill gaps in the time series [44], the spurious midpoints can be filled by the weighted average of the two nearest points. As shown in Figure 3, points a and d can be filled by the known points b and c according to Equation (3), where D ab is the time interval between points a and b, D bd is the time interval between points b and d; P a is the weight of point a; and V a is the value of point a.
The trend-fill method with a decay coefficient alpha can be used to fill the gaps in the edges of a time series. If there are fewer than five deficient points in one edge of the time series, they can be predicted according to Equation (4), where delta is the difference between points b and a and alpha is the change ratio.
three gap-fill algorithms presented in this paper. The weighted average method fill gaps in the time series [44], the spurious midpoints can be filled by the weigh age of the two nearest points. As shown in Figure 3, points a and d can be fill known points b and c according to Equation (3), where is the time interva points a and b, is the time interval between points b and d; is the weigh a; and is the value of point a. The trend-fill method with a decay coefficient alpha can be used to fill the g edges of a time series. If there are fewer than five deficient points in one edge o series, they can be predicted according to Equation (4), where delta is the diffe tween points b and a and alpha is the change ratio.
As in Figure 4, point c can be predicted from points a and b, while point predicted from points b and c. If there are more than five spurious points at th the time series, a quadratic polynomial function is used to predict the deficie Precision can be assured by using three times the number of missing points to the polynomial model. The specific gap-filling strategy follows Equation (5) As in Figure 4, point c can be predicted from points a and b, while point d can be predicted from points b and c. If there are more than five spurious points at the edge of the time series, a quadratic polynomial function is used to predict the deficient values. Precision can be assured by using three times the number of missing points to construct the polynomial model. The specific gap-filling strategy follows Equation (5), where x is the serial number of the time series point, β is the coefficient matrix, and − x is the index of the deficient point in the time series.

Enhanced SG Filtering
The SG filtering method is nearly fully automatic and only needs a few input parameters. It has been proven to be flexible and time-efficient. Its fitting index assures that the fitted time series has a low deviation from the original time series. In Equation (6), Y is the original time series value, N is the number of convoluting integers and is equal to the smoothing window size (2 m + 1), and C i is the coefficient for the ith value of the filter (smoothing window), which can be calculated from the equations presented by Madden [45]. Enhanced SG filtering methods are designed separately for time series of NDVI and other data, such as GPP and LAI. For NDVI time-series data, only smoothed points with values greater than the original value can be accepted as the upper boundary of NDVI. However, for other time-series datasets, some high values cannot be accepted and the NDVI smoothing strategy is not appropriate.
The strategy for NDVI aims to obtain the upper envelope of the time series values, as per Equation (7) [26]: where N 0 i is the original NDVI time series value after outlier recognition, N tr i is the smoothed NDVI time series value, and W i is the weight of the ith point, and d max is the maximum of the absolute difference value of N 0 i and N tr i . When the time series meets the condition: iteration times > n or abs(F K − F K−1 ) < Delta K (n is set to 20 and Delta K is set to 0.01 for the NDVI time series), then the smoothing process terminates.
In addition, for time series datasets other than NDVI, to evaluate the smoothness of the time series, the mean slope of every point in the time series is calculated. When the maximum MeanSlope of all the points meets the condition, the SG process will terminate and the specific condition can be defined as per Equation (9), where MMS is the maximum MeanSlope of all the points in the time series; MMS M is the MMS of the Mth filter process; n is the number of iterations; Delta M is the threshold for the difference between MMS M+1 and S M ; n = 5 and Delta M = 0.1.
To ensure that the smoothed time series data retains the characteristics of the original time series, all the time series are smoothed by the enhanced SG filtering method. Finally, we defined, at most, 20 smoothing times for the NDVI time series and 5 smoothing times for the GPP, reflectance, and LAI data, which can retain the original data characteristics as much as possible.

Model Fitting 2.4.1. CCTM Model
Generally, land surface changes are of four types [46]: (1) uniform variation whose variation rate is fixed; (2) decelerated variation whose absolute value of the variation rate is getting smaller and smaller; (3) accelerated variation whose absolute value of the change rate is getting smaller and smaller, which are always caused by deforestation, floods, fire, insects, urbanization, etc.; and (4) periodic variation caused by seasonality, climate variability, vegetation growth, etc. Therefore, we developed the CCTM time series model, which includes components of gradual, periodic, and abrupt variations to capture all four types of surface change. A typical model can be expressed as per Equation (10), where f 1 (y 1 ), f 2 (y 2 ), f 3 (y 3 ), f 4 (y 4 ), y 1 , y 2 , y 3 , y 4 are the functions of the Julian day, with ranges of [0, 365] or [0, 366]. In the above function-based model, f 1 (y 1 ) is used to fit the uniform variations, f 2 (y 2 ) is used to fit the accelerated variations, f 3 (y 3 ) is used to fit the decelerated variations, and f 4 (y 4 ) is used to fit the short-term periodic variations.
The typical mathematical functions can also be divided into three categories according to their second derivative ( f (x)): (1) the second derivative is larger than zero ( f (x) > 0). This kind of basic function can be a polynomial function like ax k + bx k−1 + cx k−2 + . . . + d, an exponential function like e x , or a negative logarithmic function like − ln x. (2) The second derivative is equal to zero ( f (x) = 0), such as a linear function ax + b. (3) The second derivative is lower than zero ( f (x) < 0). This kind of basic function can be an exponential function like e −x or positive logarithmic function like ln x etc.
To comprehensively represent the variations in yearly time series data, the final fitting function is composed of four parts. First, the uniform variation part can be fitted by a linear model, which can also capture the long-term change over a year: f 1 (y 1 ) = a 1,k,i y 1 + β. Second, the accelerated variation can be fitted by a positive exponential function model: f 2 (y 2 )= b e y 2 . Third, the decelerated variation can be fitted by a positive logarithmic function model: f 3 (y 3 )= c ln(y 3 ). Finally, the unimodal variation or short-term periodic variation parts can be fitted by this harmonic model [47], where f 4 (y 4 ) = d 1,k,i cos(2πy 4 ) + e 1,k,i sin(2πy 4 ). The harmonic model is comprised of a sine function and a cosine function, and it can precisely represent periodic variations. Therefore, the final model can take the form of Equation (11).
To magnify the values of the coefficients (such as a 1,k,i and b 1,k,i ), the Julian day values are normalized to the range of [0, 1] and y 4 = max−x max−min . Due to the value of exponential function being persistently larger than 0 and rapidly increasing, the Julian day has also rescaled to [0, 1] and y 2 = x−max 10 . For the logarithmic function model, the Julian days are limited to [1, e] and . Then, all of the function values are transformed to the range of [0, 1], which is beneficial for magnifying the coefficients of the fitting function.
The final CCTM time series model is as follows: where x is the Julian day; max is the maximum Julian day of the date range; min is the minimum Julian day of the date range; i is the year; a 0,i , b 0,i are the coefficients of the linear model; a 1,i , b 1,i are coefficients of the harmonic model which captures the unimodal variation; a 2,i is the coefficient of the index model that captures the accelerated variation in a year; and a 3,i is the coefficient of the logarithmic model that captures the decelerated variation in a year.

Two Patterns of the CCTM Model
Two patterns of the CCTM model have been proposed in this paper. One is the seasonal time-series pattern, another is the yearly time-series pattern. The seasonal timeseries pattern is to divide the yearly time series into four seasons, and each seasonal time series is fitted by the CCTM model, so this pattern has 24 coefficients. The yearly time-series pattern is to directly fit the yearly time series by using the CCTM model, and this pattern has 6 coefficients.
The CCTM model in the yearly time-series pattern is shown as Equation (13): In the seasonal time-series pattern, the yearly time series can be represented as the sequence of the four seasonal CCTM models (Equation (14)): where, x is the Julian day; max is the maximum day of seasonal date interval; min is the minimum day of seasonal date interval; k represents the season, where 1 = spring, 2 = summer, 3 = autumn, and 4 = winter; i is the year; a 0,i , b 0,i are the coefficients of the linear model; a 1,i , b 1,i are coefficients of the harmonic model which captures the unimodal variation; a 2,i is the coefficient of the index model that captures the accelerated variation in a year; and a 3,i is the coefficient of the logarithmic model that captures the decelerated variation in a year.

Comparison Models
To test the performance of the new model, we carried out a comparison among four models by using the above sampling points. Besides the CCTM model proposed in this paper, three other models were used, including the CCDC model, the Exp model, and the Ln model. The Exp model is the CCTM model without the decelerated variation part, and the Ln model is the CCTM model without the accelerated variation part.
The CCDC model designed for Landsat surface reflectance reconstruction provides 4, 6, and 8 coefficient function patterns, and the 6-coefficient model was applied to the comparison (Equation (15)). The Exp model is shown in Equation (16) and the Ln model is shown in Equation (17). Like the CCTM model, these models all have two patterns-a yearly time-series pattern and a seasonal time-series pattern.
where, x is the Julian date; i is the year; for all the time series datasets, the models were compared using two different preprocessing flows: (1) reconstruction without smoothing for the preprocess; and (2) reconstruction with smoothing for the preprocess.

Evaluation of Time Series Reconstruction Accuracy
The gap-filled time series has maintained the original characteristics in the time series variations, so it can be used as the standard time series for evaluating the reconstruction effect. Two indices between the gap-filled time series and reconstructed time series are calculated to evaluate the reconstruction effects. First, the RMSE (Root mean square error index) can measure the degree of deviation of the model from the original data. Second, the R 2 index can evaluate the consistency of the original data before and after reconstruction. These indices are shown as Equation (18)

Accuracy Validation of Yearly and Seasonal Time Series Patterns
The continuous reconstruction results of the four models are compared in the two patterns: (1) yearly time series reconstruction; and (2) seasonal segment time series reconstruction. Both continuous-time series reconstructions use the above-mentioned methods including the time series outlier gap-filling algorithm and the Enhanced SG filter noise smoothing process, and then the accuracies of the four models for the continuous-time series reconstruction are evaluated by R 2 and RMSE respectively.
For the yearly time series pattern, the following steps compare the reconstruction results of the four models in the yearly time series reconstruction pattern. First, the yearly cropland, forest, and grassland NDVI time-series datasets are selected to evaluate the trend-fitting ability for its application in phenology extraction. In total, 494 cropland sampling points, 456 forest sampling points, and 2261 grass sampling points are chosen to evaluate the fitting abilities among different models. All the NDVI time series points are preprocessed, including outlier recognition, Gap filling, and the Enhanced SG filter. The yearly time series that remove some detailed variations and only retains the yearly trend is not too challenging for the four models. Therefore, for the yearly time series, the SG filter iteration time is set to be less than 100, and the yearly characteristics will not be lost but detailed variations are nearly lost during the process. Six or 8 parameters are not too challenging for smoothed time series. Second, the fitting accuracies are evaluated using the indices R 2 and the RMSE calculated from the preprocessed time series and the fitted time series.
For the seasonal time series pattern, the time series of the sampling points in Section 2.3.1 are preprocessed by outlier recognition and gap-filling methods in Sections 2.3.1 and 2.3.2. The preprocessed time series was then smoothed by the Enhance SG filter. Finally, with the seasonal division strategy, the time series are separately reconstructed by the four function fitting models. The fitting R 2 and RMSE are separately calculated from the preprocessed time series and the final time series to evaluate the fitting accuracy of the four function fitting models.

Comparison of Yearly Timeseries Pattern
The RMSE probability distribution plots in Figure 5 and the R 2 probability distribution plots in Figure 6 show that the Harmonic model has the strongest fitting ability for the yearly time series. More than 95%R 2 values and RMSE values are separately larger than 0.90 and less than 0.01 for the Harmonic model and the CCTM model. The CCTM model is slightly inferior to the Harmonic model in fitting accuracy, and the Exp model and the Ln models have shown nearly the same fitting ability. The average statistical result of R 2 and RMSE has shown that the Harmonic model has the best fitting accuracy, which has minimal differences from the original time series.
In addition, Figure 7 shows some of the fitting curves of different models in yearly time series patterns. The four models have some differences in the fitting effects. Although the Harmonic model and the CCTM model have nearly the same fitting effects, the Harmonic model has shown more time-series symmetrical characteristics than the CCTM model. The annual time series contains more periodic variation characteristics than the seasonal time series, and to comprehensively assess the fitting effects of the four models, the seasonal segment time-series pattern will also be compared among the four function fitting models. the seasonal segment time-series pattern will also be compared among the four function fitting models.   the seasonal segment time-series pattern will also be compared among the four function fitting models.   the seasonal segment time-series pattern will also be compared among the four function fitting models.

Smoothing Process
Before the smoothing process of the SG filter, the gap-filled time series are extracted, and R 2 and RMSE between the gap-filled time series and smoothing processed time series are calculated. The final statistical results are shown in Table 2. For the NDVI data, the smoothing process was used to reduce noise. Nearly all of the data except for shrubland have RMSEs > 0.01, and almost 95% of data except that of EBF have R 2 -values >0.9.
For reflectance data, more than 90% of time series points have RMSEs <0.035, and more than 90% of land cover types except shrubland have R 2 -values >0.9. The data for forest types, including the ENF, DBF, and MF, have higher R 2 -values than other land covers. In addition, the reflectance data of land cover types except shrubland have lower RMSEs and higher R 2 -values. This can, to some extent, indicate that the smoothing process reduces noise in the time series data and maintains its original long-term trend.

ENF
For the LAI dataset (Table 3), the harmonic model shows the best reconstruction performance, with the CCTM model being slightly inferior. Except for the EBF land cover type, nearly 90% of time series points have fitting R²-values >0.9, which indicates that the CCTM model also has a great trend-preserving ability. For different land cover types, the For the GPP dataset (Table 3), more than 95% of the time series points of land cover types except EBF have fitting R 2 -values >0.98, and more than 95% of time series points of all the land cover types have fitting RMSEs <0.02 g C m −2 using the CCTM model. The CCTM model is significantly better than other models, especially for the GPP time series of shrubland, savanna, grassland, cropland, and EBF. For other forest types, except for EBF, there are smaller differences among the four models. Overall, the CCTM model performs well for the GPP time series reconstruction for all land cover types.
For the LAI dataset (Table 3), the harmonic model shows the best reconstruction performance, with the CCTM model being slightly inferior. Except for the EBF land cover type, nearly 90% of time series points have fitting R 2 -values >0.9, which indicates that the CCTM model also has a great trend-preserving ability. For different land cover types, the reconstruction RMSEs are relatively inconsistent, and the forest types have larger RMSEs than non-forest types such as savanna, shrubland, and cropland. Due to the relatively low RMSE and high R 2 -values, the CCTM model is suitable for LAI time series reconstruction.
For the reflectance dataset, the reconstruction result after smoothing achieved high R 2 -values (nearly 100% of time series points have R 2 -values >0.9 by any of the models)( Figure 10). The CCTM model has the best reconstruction effect among the four models, and nearly all the time series points have R 2 -values >0.95( Figure 11). In addition, the CCTM also gained the lowest RMSE, about 90% of time series points of all land cover types have RMSE <0.01.

Discussion
The discussion section consists of three parts: (1) the major advantages and limitations, the reasons for high accuracy reconstruction results and the limitations of the new method are discussed in this part; (2) the potential applications; and (3) the time series functionalization and compression, the compression patterns and choices are discussed in this section.

Discussion
The discussion section consists of three parts: (1) the major advantages and limitations, the reasons for high accuracy reconstruction results and the limitations of the new method are discussed in this part; (2) the potential applications; and (3) the time series functionalization and compression, the compression patterns and choices are discussed in Figure 11. R 2 probability distribution of reconstructed MSR time series data after smoothing for different landcover types. The frequency is the number of time series points in a certain RMSE range, and different landcover types may have different RMSE categories.

Discussion
The discussion section consists of three parts: (1) the major advantages and limitations, the reasons for high accuracy reconstruction results and the limitations of the new method are discussed in this part; (2) the potential applications; and (3) the time series functionalization and compression, the compression patterns and choices are discussed in this section.

Advantages and Limitations
In this paper, the fitting abilities of the four function fitting models are compared in the two time series patterns-the yearly time series pattern and the seasonal time series pattern. Different results exist in the two different patterns. For the yearly time series pattern, the Harmonic model shows the best fitting accuracy, and the CCTM model is slightly inferior to the Harmonic model. Most of the cropland, forest, and grassland NDVI time series have obvious seasonal time-series variation characteristics [48][49][50], and the Harmonic model has the most symmetrical fitting functions among the four function fitting models [51], so the Harmonic model has shown the best fitting accuracy among the four models.
For the seasonal time series pattern, the four function fitting models have all shown higher fitting accuracy than the models in the yearly time series pattern. Obviously, the seasonal pattern is a kind of segment fitting technique, and therefore the higher fitting accuracies were obtained in this pattern. The CCTM model has shown the highest fitting accuracy in this pattern. The CCTM model is designed according to the type of time series variation rates. It is a more comprehensive function fitting model than the Harmonic model. When it comes to time series with fewer periodic variations, such as the time series in one season, the CCTM model shows better fitting ability than the Harmonic model.
The higher accuracy of the CCTM model also shows that the CCTM model has a stronger comprehensive function fitting ability than the Harmonic model for the asymmetric time-series variations induced by vegetation growth, natural disasters, etc. From the meaning of the second derivative mathematical function, the CCTM model has more comprehensive representation techniques. Therefore, the CCTM model can be more accurate than the Harmonic model for the continuous-time series reconstruction.
The seasonal segment fitting results show that there are differences among different vegetation types. Obviously, different vegetation types mean different seasonal variation cycles and ranges of time series values, which lead to different fitting effects among different function fitting models. For the EBF vegetation type, the periodic variation characteristic is not significant in the yearly NDVI time series, and the Harmonic model and the CCTM model showed the largest differences.
However, the CCTM model has more complex function fitting structures than the Harmonic model. If there is no requirement for high-precision reconstruction results or model interpretability, the CCTM model is not suggested.

Potential Applications
Due to the high precision reconstruction effect, the CCTM model may also be applied to land surface phenology detection and monitoring, seasonal characteristics extraction, and trend analysis. For seasonal characteristic extraction, the CCTM model is specially designed to decompose the seasonal characteristics into four change parts and effectively extract the seasonal variation characteristics. For the trend analysis, because the CCTM model can decompose the time series into different variation types, the CCTM model can analyze trends of different variation types. The phenology extraction demands extracting the time series statistical features to obtain different transition states [52]. The CCTM model can represent the time series practically with high accuracy, improving the efficiency and accuracy of phenological extraction.

Time Series Functionalization and Compression
The CCTM model also provides a solution to minimize time series storage, one-year time series can be precisely characterized with only 24 coefficients in the seasonal time series pattern. The 24 coefficients can replace the original huge dataset as the new input, conducive to analyzing and applying extensive time-series data. Because extracting statistical features with a function is easier than extracting from actual time series, the functionalization for the time series can also optimize the extraction efficiency of statistical characteristics.
In addition, the coefficients have mathematical meanings, and the coefficients of different variation types represent the weights of different variations. More research is needed to explore the correlation of the coefficients with different variation types.
The CCTM model can also be applied to time series image compression. The daily time series have 365 or 366 bands, the 4-day time series dataset has approximately 90 bands, and the 8-day time series dataset has about 46 bands. All these time series can be stored in only 24 bands in the seasonal time series pattern, with each seasonal time series being fitted by the CCTM model which has 6 coefficients. Therefore the 4-day time series can approximately reach the compression ratio of 3.75, the 8-day time series dataset can approximately achieve the compression ratio of 1.91, and the daily time series dataset can approximately achieve the compression ratio of 15.
The DBF land cover type's daily reflectance dataset and 8-day NDVI dataset have been selected to evaluate the compression effects. The CCTM model separately compresses these two datasets with 24 coefficients in seasonal time series pattern and 6 coefficients in yearly time series pattern.
The compression of the CCTM model has two patterns. The seasonal time series compression pattern which has 24 coefficients is to divide the yearly time series into four seasons. The yearly time series compression pattern which has 6 coefficients is to directly fit the yearly time series by using the CCTM model. By using the above compression patterns, the storage before or after compression is compared and the R 2 between the original time series and compressed time series is calculated ( Table 4). The original NDVI compression image has 383 rows, 307 columns, and 46 bands. Every pixel value is stored with a 32-bit integer, so the total image size is 20.63 MB. Similarly, the original consistent reflectance image has 365 bands, so the total image size is 163.7 MB. For the seasonal compression pattern, when applying the CCTM model to these datasets, the storage of the compressed time-series dataset for both NDVI and reflectance time series is only 10.7 MB, which is 52% of the original 8-day NDVI dataset and 6.5% of the original daily reflectance dataset. For the yearly time series compression pattern, the final compressed NDVI and reflectance time series can be stored with only 2.69 MB, which is 13% of the original 8-day NDVI dataset and 1.6% of the original daily reflectance dataset. The number of compression coefficients should be chosen using the actual data precision when selecting the compression method.
The CCTM model is flexible on the number of the coefficients; a whole year time series can be compressed with 6 and 24 parameters. Remote sensing image compression has its standards, and the specified compression methods using the CCTM model need further research.
Due to the seasonal segment time series fitting technique, the CCTM model can achieve higher accuracy in the seasonal time series pattern than that in the yearly time series pattern, but the compression in the seasonal time series pattern needs more storage space than the compression in the yearly time series pattern. The compression pattern should be decided by the trade-off between compression ratio and compression accuracy. Only when the yearly time series compression pattern can not meet the accuracy requirement, the seasonal time series compression pattern should be adopted.

Conclusions
The essence of remote sensing time series reconstruction is to repair missing values by exploring and using relevant time series data. Traditional time series reconstruction methods commonly have over-smoothing and under-fitting deficiencies, and traditional reconstruction function-based models are too simple to reflect complex time series variations. Certainly, a complex model with many more parameters definitely fits data better than a simple model, but simply increasing the number of parameters is not effective. The CCTM model with 6 parameters is constructed based on the second-order derivative theorem and this is why its reconstruction precision is higher than the Harmonic model with 8 parameters. In general, the CCTM model can comprehensively reflect the complex characteristics of time series variation and embody seasonal characteristics.
Compared to the Exp, Ln, and harmonic models, the CCTM model has lower fitting errors. In the reconstruction of GPP and MSR time series with a smoothing process, >95% of time series points (except the EBF type) had R 2 -values >0.9 and RMSEs < 2 g C m −2 and 0.025, respectively by using the CCTM model. For NDVI reconstruction after smoothing, >95% of time series points had R 2 -values >0.9 and RMSEs < 0.01 (except for the ENF and MF types) using the CCTM model. The CCTM model had the fewest extreme fitting errors but also has its limitations. Because the TP and LST time-series datasets contain too many frequent fluctuations, none of the models is suitable for their reconstruction.
Developing high-precision models can improve the precision of time series reconstruction and be beneficial in obtaining accurate knowledge of the variations in the time series dataset. Besides, the CCTM model is supposed to be applied to time series change detection, trend analysis, extraction of phenology and seasonal characteristics, time-series image compression, and land cover classification.