*3.2. Adaptive Smoothing*

The proposed reconstruction algorithm requires complete time series with contiguous observations in the temporal domain for each spectral band. Gaps due to missing or masked observations are filled with linearly interpolated values of clear observations. In case values are masked near the beginning and end of the time series, the nearest clear observation is extrapolated.

Even after interpolating the masked values, the reflectance values R*λ*,*<sup>i</sup>* (with *λ* representing the wavelength and *i* the acquisition time) can still be noisy. For instance, the acquisition of 7 June 2019 corresponds to a "clear" observation (vertical line 2 in Figure 3). However, the reflectance in bands B4 and B8 is higher than expected, based on the previous and next clear observations (see Figure 3d). The isolated high reflectance can be interpreted as an outlier. This is also confirmed by the corresponding NDVI value, N*<sup>i</sup>* (see Figure 3a). Similar to the time series of the interpolated spectral reflectance, R*λ*,*i*, the time series for the interpolated NDVI value, N*i*, also shows an outlier for the acquisition from 7 June 2019. The isolated drop in NDVI is not compatible with a gradual change of vegetation and can be regarded as noisy. By inspecting Figure 3a, it is shown that the image scene indeed contains some cloud cover that could explain the noise.

Several algorithms for smoothing time series have been proposed in the literature, including the Whittaker smoother [40], the Savitzky–Golay filter [41], and its modifications [20,42]. Because the proposed RTSR algorithm was inspired by the work in [20], the same smoothing filter (Savitzky–Golay) was chosen here. It is important to note that other smoothing algorithms can easily be adopted in the proposed reconstruction method. Similar to the approach in [20] to reconstruct the NDVI time series, a first Savitzky–Golay filter, referred to as SG1, creates the long-term change trend (represented by the orange lines in Figure 3). It has a relatively wide half-width of the smoothing window (7). This value is based on our own experiments and on values found in the literature [20]. Likewise, the parameter for the polynomial degree was set to 2 for SG1 (suggested between 2 and 4 in [20]). Larger values are supposed to better follow higher frequencies in the time series.

As illustrated in Figure 3b, the long-term change trend of the NDVI, SG1(N*i*), underestimates the NDVI values for the clear observations in May. Based on the assumption that relatively high NDVI values correspond to trustworthy observations, the upper NDVI envelope is expected to better reflect the dynamic changes of interest in the NDVI temporal profile [20–22]. The acquisitions with a calculated NDVI above the long-term change trend (green dots) are therefore assumed to be more reliable than those with a calculated NDVI value below the long-term change trend (orange dots). Following this reasoning, the original spectral reflectance values R*λ*,*<sup>i</sup>* observed in May can be considered as reliable. Unlike the NDVI, however, the reflectance values are lower than the long-term change trend SG1(R*λ*,*i*) (indicated in orange in Figure 3e). The assumption that clouds or poor atmospheric conditions depress NDVI values does not hold for the individual spectral bands. For instance, reflectance values in red and near infrared typically increase in case of cloud cover, but decrease in case of cloud shadow. Neither the upper nor the lower envelope reflect the dynamic changes of interest in the spectral reflectance temporal profile that would be observed in clear conditions. A different algorithmic approach for spectral reflectance was therefore developed, which is listed in Algorithm 1.

The reconstructed reflectance time series Rˆ *<sup>λ</sup>*,*<sup>i</sup>* is initialized as the long-term change trend SG1(R*λ*,*i*). If, for a given acquisition *i*, the interpolated NDVI value N*<sup>i</sup>* (calculated from the interpolated reflectance R*λ*,*i*) shows a higher value than the current NDVI estimate Nˆ *<sup>i</sup>* (calculated from the estimated reflectance Rˆ *<sup>λ</sup>*,*i*), then the corresponding interpolated reflectance value R*λ*,*<sup>i</sup>* is assumed to be reliable and replaces the current estimate Rˆ *<sup>λ</sup>*,*<sup>i</sup>* (green dots in Figure 3e). Else, the current estimate Rˆ *<sup>λ</sup>*,*<sup>i</sup>* is retained in favor of the interpolated reflectance value R*λ*,*<sup>i</sup>* that is assumed to be noisy for the given acquisition *i* (e.g., orange dot in Figure 3e). By comparing the green and orange dots in both Figure 3b,e, it is shown that the decision which value to select (smoothed or interpolated version) is driven by the NDVI time series.

The algorithm proceeds by smoothing the retained spectral values with a second version of the SG filter, referred to as SG2. The smoothed values are then used for the next iteration. The reconstructed spectral reflectance time series Rˆ *<sup>λ</sup>*,*<sup>i</sup>* after the final iteration is shown in Figure 3c,f in blue. Different stopping criteria for number of iterations can be implemented. In [20], a fitting-effect index was proposed. Here, the number of iterations was experimentally fixed to five. More iterations did not further improve results. As a comparison, the reconstructed time series based on the dynamic temporal smoothing (DTS [33]) method is also shown (in green).

**Algorithm 1** Reflectance reconstruction algorithm Rˆ *<sup>λ</sup>*,*<sup>i</sup>* <sup>←</sup> SG1(R*λ*,*i*) initialize reconstructed reflectance with long-term change trend **while** Stopping criterion not met (e.g., 5 iterations) **do if** Nˆ *<sup>i</sup>* < N*<sup>i</sup>* **then** observation is clear Rˆ *<sup>λ</sup>*,*<sup>i</sup>* <sup>←</sup> <sup>R</sup>*λ*,*<sup>i</sup>* update reconstructed reflectance **end if** Rˆ *<sup>λ</sup>*,*<sup>i</sup>* <sup>←</sup> SG2(Rˆ *<sup>λ</sup>*,*i*) smoothen reconstructed reflectance **end while**
