## 2 Materials and Methods
### Phenological Metrics Extraction  
Satellite remote sensing time series are affected by environmental conditions such as clouds and atmospheric contaminants, and the viewing angle of the satellite instruments can introduce additional uncertainties. To enhance the seasonal variations of interest indices, it is necessary to apply cleaning and curve-fitting methods to enhance the annual variations of the time series. Before fitting, data cleaning and smoothing were performed to reduce the noise influences. First, the obtained time series were filtered to remove outliers using the method proposed by Timesat 3.3 [37]. When iterating over the time series, the following criteria were used to remove outliers: data points that exceeded 2 standard deviations (named as a cutoff) away from the mean of the time series, less than the average of its nearest validate values in the time series (named as immediate neighbors) minus the cutoff, or higher than the maximum of its immediate neighbors plus the cutoff. Fourier transform (FT), a basic signal-processing method, is a traditional signal noise reduction technique that conducts frequency domain analysis and has been applied in a number of phenological studies [38; 39]. Evenly spaced time series are required in FT; therefore, the pre-processed time series were binned into 10-day windows and the average value was calculated to characterize the signal. Linear interpolation of the immediate neighbor windows was used to fill gaps in the time series. If satellite data was consecutively missing for more than 30 days (which occurred in the S2 data of 2017 and 2018 in this research), the data before the missing period would be discarded to reduce the effect of linear interpolation on the curve-fitting accuracy. The FT decomposes signals into a series of harmonics sorted by frequency and one base value, in which non-periodic noise is usually characterized by high-frequency harmonics and seasonal trends are usually characterized by harmonics of lower frequency [40]. The non-periodic noise can be removed after inverse FT is conducted using the first few low-frequency harmonics and the base value [41]. [PERSON] et al. [42] reported that with the first four harmonics and the base value, the phenological profile can be reconstructed reasonably, and this study conducted the inverse FT in the same way. The DL function was constructed using two logistic functions to simulate the single peak phenology. Its basis function can be written as [43]:  
\[\mathrm{g}(\mathrm{t},a_{1},a_{2},b_{1},b_{2},\mathrm{c},\mathrm{d})=\mathrm{ c}+\mathrm{d}\bigg{(}\frac{1}{1+e^{b_{1}(a_{1}-t)}}-\frac{1}{1+e^{b_{2}(a_{2}-t)}} \bigg{)} \tag{3}\]where \(a_{1},a_{2}\) are the inflection points of the growing period and falling period, respectively, \(b_{1},b_{2}\) are the slopes of \(a_{1},a_{2}\), respectively, c is the minimum value of the curve, and d is the amplitude of the curve. A non-linear least square function was considered in the curve fitting. Each DL function can fit the time series of polarization combinations or NDVI for one growing season, approximately one year for forests, so the time series were divided by year during the DL curve fitting. The FT, the inverse FT, and the curve fitting process of the DL were performed in Python 3.8. The FT and the inverse FT were conducted by NumPy [44], while the curve fitting was conducted by a least-squared optimization method from SciPy [45]. To facilitate the comparison of phenological phases extracted by this research and the phenological metrics in the PhenoCCam Level 5 product, the same dynamic threshold method was applied to extract the phenological metrics, for each growing period or falling period, 50% of its amplitude was set as the threshold [26]. Figure 2 shows the curve fitting and the extracting of phenological metrics. The first day, when the data increase exceeded the threshold, was marked as SOS, and the first day when the data decreased below the threshold was marked as EOS.