Variations in Phenology Identification Strategies across the Mongolian Plateau Using Multiple Data Sources and Methods

Abstract

Satellite data and algorithms directly affect the accuracy of phenological estimation; therefore, it is necessary to compare and verify existing phenological models to identify the optimal combination of data and algorithms across the Mongolian Plateau (MP). This study used five phenology fitting algorithms—double logistic (DL) and polynomial fitting (Poly) combined with the dynamic threshold method at thresholds of 35% and 50% (DL-G35, DL-G50, Poly-G35, and Poly-G50) and DL combined with the cumulative curvature extreme value method (DL-CUM)—and two data types—the enhanced vegetation index (EVI) and solar-induced chlorophyll fluorescence (SIF)—to identify the start (SOS), peak (POS), and end (EOS) of the growing season in alpine meadow (ALM), desert steppe (DRS), forest vegetation (FV), meadow grassland (MEG), and typical grassland (TYG) of the MP. The optimal methods for identifying the SOS, POS, and EOS of typical grassland areas were Poly-G50 (NSE = 0.12, Pbias = 0.22%), DL-G35/50 (NSE = −0.01, Pbias = −0.06%), and Poly-G35 (NSE = 0.02, Pbias = 0.08%), respectively, based on SIF data. The best methods for identifying the SOS, POS, and EOS of desert steppe areas were Poly-G35 (NSE = −0.27, Pbias = −1.49%), Poly-G35/50 (NSE = −0.58, Pbias = −1.39%), and Poly-G35 (NSE = 0.29, Pbias = −0.61%), respectively, based on EVI data. The data source explained most of the differences in phenological estimates. The accuracy of polynomial fitting was significantly greater than that of the DL method, while all methods were better at identifying SOS and POS than they were at identifying EOS. Our findings can help to facilitate the establishment of a phenological estimation system suitable for the Mongolian Plateau and improve the observation methods of vegetation phenology.

1. Introduction

Vegetation phenology refers to the different timing nodes (e.g., germination, leaf unfolding, flowering, fruiting, senescence, and dormancy) that occur periodically in vegetation over many years [1,2]. Vegetation phenology is an important indicator of climate change and plays a vital role in regulating ecosystem structures and functions [3]. Land surface phenology (LSP) observes such seasonal changes in land surface vegetation through sensors [4] and contains three commonly used, important indicators, namely, the start (SOS), peak (POS), and end (EOS) of the growing season. Previous studies on LSP have been limited to comparing phenology data sources or algorithms, and differences in phenology identification caused by diverse combinations of data sources and methods remain unclear. Phenology plays an important role in maintaining the stability of the overall carbon–water cycle of the ecosystem [5]; therefore, selecting appropriate phenology identification schemes for different ecosystems and accurately estimating large-scale phenology play key roles in ecological and environmental management.

Satellite remote sensing has been extensively used in large-scale phenological research, and many studies have shown significant differences in phenological estimates based on different remote sensing data sources [6,7]. Current phenology estimation data sources are primarily divided into two categories: (1) The greenness of vegetation canopy morphology, represented by the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) [8]. Obtaining vegetation phenology information from the perspective of vegetation morphology and greenness using such indices has been widely adopted [9]. However, some studies have found that spring starts earlier, leaf shedding occurs later [10,11], and the growing season is longer when identification is based on VI data [12]. (2) The vegetation physiological signal, represented by solar-induced chlorophyll fluorescence (SIF). As it does not contain nonfluorescent targets, SIF has unique advantages in characterizing the phenological information of vegetative physiological development [13]. SIF has been increasingly applied in phenological studies in areas with sparse vegetation coverage [14,15]. However, comparisons of phenological identification based on these two data sources and evaluations of their applicability in different ecosystems are lacking.

Climate change has led to the degradation of ground vegetation in some areas, creating new challenges in the identification accuracy of current phenological algorithms [16]. Phenological studies have attempted to determine phenological estimation models suitable for arid regions with sparse vegetation [17]. The present general LSP extraction process is primarily divided into two steps. The Savitzky–Golay filter [18], asymmetric Gaussian filter [19], double logistic function [20], time-series harmonic analysis [21], and other methods have been used for the first step: reconstructing a time-series curve. The threshold [22], derivative [23], fitting [24], and other methods have been used for the second step: extracting key phenological parameters. Previous research has focused on the comparison of two aspects: phenological parameter extraction algorithms [25] and phenological parameters calculated from diverse data collections [26,27]. However, limited attention has been paid to the aspects and attributes of vegetation phenology that reflect the choice of methods and data. Therefore, current phenological research aims to compare the accuracy of phenological identification based on different data and methods and determine the optimal combination of data and algorithms for different vegetation types.

The Mongolian Plateau (MP) is located in an arid/semiarid region with a vast grassland area and is an important ecological barrier in northern Asia. In the context of severe global climate fluctuations, changes in plant phenology will affect the energy exchange in the Eurasian ecosystem and feedback to the global climate system [28]. However, most phenological identification studies considered differences in phenological estimation in isolation from data [29] or methods [30], but the primary factors affecting the accuracy of phenological estimations could not be identified. It is impossible to objectively evaluate phenological characteristics on a large scale. Therefore, this study analyzed the differences in the temporal and spatial characteristics of phenology on the MP by identifying three phenological parameters using different data sources and methods.

2. Materials and Methods

2.1. Study Area

The main part of the MP, including Mongolia and the Inner Mongolia Autonomous Region of northern China, was selected as the study area (Figure 1, 87°40′E–122°15′E, 37°46′N–53°08′N). This region has apparent seasonal variations in natural vegetation, a fragile and sensitive ecological environment, complex and variable physical geography, and diverse landscape ecological types. The elevation increases gradually from east to west, with an average of 1580 m (Figure 1a). The vegetation cover types include alpine meadow (ALM), desert steppe (DRS), forest vegetation (FV), meadow grassland (MEG), and typical grassland (TYG) (Figure 1b).

2.2. Data Sources

2.2.1. Eddy Covariance Flux Data

Daily gross primary productivity (GPP) from flux data based on the eddy covariance (EC) method of the China Flux Observation Consortium (ChinaFLUX, http://www.chinaflux.org/, accessed on 11 May 2023) were compared with LSP estimated from remote sensing data (Figure 2). Three sites located in the study area (Figure 1b) were selected for analysis (

Table 1. Flux site names, coordinates, and land cover type from 2003 to 2018.

Site Name	Lat. (°N)	Lon. (°E)	Land Cover Type	% of Land Cover	Available Range
Baiyin Xile	43.326	116.404	Typical grassland	21.1%	2004–2010
Duolun	42.047	116.284	Typical grassland	21.1%	2006–2010
Siziwang	41.790	111.897	Desert Steppe	28.4%	2010–2012

% of total land cover = area of one land cover type/total area of the study area*100%.

): two TYG sites (Baiyin Xile and Duolun, Inner Mongolia) and one DRS site (Siziwang, Inner Mongolia).

We collected daily GPP data from the study field sites. The data products were calibrated using the method described by Papale et al. [31], followed by gap filling and partitioning as described by Reichstein et al. [32]. Net ecosystem exchange (NEE) was divided into ecosystem respiration (Re) and GPP after quality control, where GPP was calculated using Equation (1):

$$\mathrm{G}\mathrm{P}\mathrm{P}=\mathrm{R}\mathrm{E}-\mathrm{N}\mathrm{E}\mathrm{E}$$

(1)

2.2.2. Remote Sensing Data

EVI data from 2001 to 2020 were derived from the 16-day MOD13A1 (Collection 6.1) EVI product provided by NASA (https://e4ftl01.cr.usgs.gov/, accessed on 20 October 2022) with a spatial resolution of 500 m; this data type is widely used to infer vegetation phenology and growth. The data are calibrated for errors caused by adverse atmospheric, radiation, and geometric conditions and can be used to characterize the vegetation radiation and greenness status. In this study, the spatial resolution was resampled to 0.05° to match the fluorescence dataset. The MOD13A1 EVI data have large areas with filling values or pixels of poor quality (such as cloud-coverage pixels), which lead to a large amount of missing data and affect experimental results. When performing quality control on the time series, cloud-coverage pixels (to be reconstructed) were filled with the weighted average of adjacent pixels of better quality (reference pixels). The weight was determined by the spectral angle similarity (distance between cells) between the pixel to be reconstructed and the reference pixel [33]. Finally, the amplitude-limit method was used to greatly increase the proportion of pixels that can obtain phenology in the MP [17]. The “amplitude-limit method” refers to a technique used to restrict the magnitude or range of certain values in a given context and is utilized in situations where extreme or unrealistic values must be constrained to ensure the validity and reliability of results. That is, when the number of pixels of >0.09 in the specified time series was >4, phenology was extracted; otherwise, it was eliminated.

SIF data were obtained from the GOSIF dataset (https://globalecology.unh.edu, accessed on 21 June 2022) based on the inversion of the OCO-2 global SIF product provided by Xiao et al., with a high spatiotemporal resolution of continuous global coverage and long-term records (2001–2020, 8-day, 0.05°). SIF effectively characterizes the photosynthetic state of vegetation [34]. After quality control, according to the fluctuation characteristics of the curve, pixels whose peak value in the growing season was >60% of the maximum value before turning green were used for phenological extraction. Otherwise, they were excluded from the analysis.

Land-cover-type data were obtained from the MCD12C1 land cover product (https://lpdaac.usgs.gov/products/mcd12c1v006/, accessed on 20 October 2022) with a spatial resolution of 0.05° and were used to analyze differences in phenology retrieval algorithms for different vegetation types.

2.3. Methods

2.3.1. Phenology Extraction Algorithms

Five methods were applied as phenology extraction algorithms to SIF and EVI data to identify the SOS, EOS, and POS from 2001 to 2020. This study used five phenology fitting algorithms—double logistic (DL) and polynomial fitting (Poly) combined with the dynamic threshold method at thresholds of 35% and 50% (DL-G35, DL-G50, Poly-G35, Poly-G50) and DL combined with the cumulative curvature extreme value method (DL-CUM). The unit was recorded as the day of the year (DOY), and the combinations of algorithms were as follows:

(1)

DL-G35

This method was based on fitting time-series curves using dual logistic fitting functions and then using dynamic thresholds defined as 35% to determine phenological parameters (Figure 3a). SOS was defined as the first day when the EVI_ratio value exceeded 0.35 (G35) in a year, and EOS was defined as the first day after SOS when the EVI_ratio value fell below 0.35 (G35) in a year. POS was the date of the peak of the fitted curve [22]. The specific algorithms were as follows:

$$\mathrm{y}\left(\mathrm{t}\right)=\frac{\mathrm{c}}{1+{\mathrm{e}}^{\mathrm{a}+\mathrm{b}\mathrm{t}}}+\mathrm{d}$$

(2)

$${\mathrm{E}\mathrm{V}\mathrm{I}}_{\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}}=({\mathrm{E}\mathrm{V}\mathrm{I}}_{\mathrm{t}}-{\mathrm{E}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}})/({\mathrm{E}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{a}\mathrm{x}}-{\mathrm{E}\mathrm{V}\mathrm{I}}_{\mathrm{m}\mathrm{i}\mathrm{n}})$$

(3)

where EVI_t is the EVI value for a given time t and EVI_max and EVI_min are the maximum and minimum EVI values in the annual EVI cycle, respectively.

(2)

DL-G50

This method was the same as DL-G35 except that the dynamic thresholds were defined as 50% and the corresponding critical EVI_ratio value was 0.5 (G50).

(3)

Poly-G35

Time-series curves were fitted using a polynomial fitting function (Equation (4)), and phenological parameters were extracted using the threshold values shown in Figure 2 [35]:

$$y=a_6{x}^6+a_5{x}^5-a_4{x}^4+a_3{x}^3-a_2{x}^2+a_{1}x+a_0$$

(4)

where x is the DOY and a is the fitting coefficient. A binary six-polynomial fit was used.

(4)

Poly-G50

This method was the same as Poly-G35 but used the 50% threshold.

(5)

DL-CUM

This method extracts vegetation phenological parameters based on the cumulative vegetation index [36]. The original SIF/EVI time-series curve was converted into a “prologue-type NDVI cumulative curve.” The DL function was then used to fit the cumulative NDVI time curve, and the phenological parameters were extracted according to the curvature extremum K. As shown in Figure 3b, SOS was defined as the time when K reached the first derivative maximum, EOS was the time of the first derivative minimum in the late growing season, and POS was the point at which the positive and negative values of K changed. The formulae used were as follows:

$$\mathrm{K}=\frac{\mathrm{d}\mathsf{\alpha }}{\mathrm{d}\mathrm{s}}=-\frac{{\mathrm{b}}^2\mathrm{c}\mathrm{z}(1-\mathrm{z}){(1+\mathrm{z})}^3}{{[{\left(1-\mathrm{z}\right)}^4+{\left(\mathrm{b}\mathrm{c}\mathrm{z}\right)}^2]}^{1.5}}$$

(5)

$$\mathrm{z}={\mathrm{e}}^{(\mathrm{a}+\mathrm{b}\mathrm{t})}$$

(6)

where t is the date; y(t) is the EVI at date t; a and b are the estimated empirical coefficients; [c + d] is the maximum value cumulative EVI or SIF; and d is the minimum value, i.e., background value.

2.3.2. Evaluation of Land Surface Phenology Retrieval Algorithm

The phenology indicators were extracted from the daily GPP time-series data of each ChinaFLUX site using the phenology identification algorithm corresponding to the remote sensing data and were used as a standard for judging the accuracy of phenology identification. These phenology indicators were compared with the LSP indicators generated using the EVI and SIF data. Nash–Sutcliffe efficiency (NSE) [37] and percentage bias (Pbias) [38] were used as indicators of accuracy when assessing the agreement between the LSPs generated from the remote sensing data (EVI and SIF) and the GPP data. R² was not used because of large uncertainties in its application to event-based simulations with small sample sizes [39]. Calculation of the indices used in this study was as follows:

$$\mathrm{N}\mathrm{S}\mathrm{E}=1-\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{s},\mathrm{i}}-{\mathrm{O}}_{\mathrm{m},\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{m},\mathrm{i}}-\overline{{\mathrm{O}}_{\mathrm{m}}}\right)}^2}$$

(7)

$$\mathrm{P}\mathrm{b}\mathrm{i}\mathrm{a}\mathrm{s}=\frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{O}}_{\mathrm{s},\mathrm{i}}-{\mathrm{O}}_{\mathrm{m},\mathrm{i}}\right)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{O}}_{\mathrm{m},\mathrm{i}}}$$

(8)

where n is the number of simulations; O_s,i is the simulated values output by the algorithms based on SIF and EVI data in the i-th year; and O_m,i is the average output value of the measured data by GPP in the i-th year. NSE is in the range from −∞ to 1, and the closer the value is to 1, the closer the simulation result is to the observed data, and the better the model performance. A Pbias with a low magnitude indicates an accurate model simulation, and a Pbias of zero indicates no systematic deviation between the simulated and observed values. Positive values indicate model underestimation bias and negative values indicate model overestimation bias.

3. Results

3.1. Evaluation of Land Surface Phenology Retrieval Data and Algorithms

The identification results of phenological parameters based on EVI data were verified using the site data of three flux tower sites for a total of 15 years (Figure 4). The results showed that the simulated phenological parameters in the typical steppe were closer to the phenological results of the flux GPP than those in DRS. Specifically, for TYG, the Poly-G50 (NSE = −0.53, Pbias = −0.57%) method gave the most accurate estimate of SOS (Figure 4d); DL-CUM significantly overestimated SOS and DL-35 significantly underestimated it. Estimating POS by DL fitting was optimal, and DL-CUM had the smallest deviation (NSE = 0.04, Pbias = 0.19%). EOS from DL-CUM (NSE = 0.28, Pbias = −0.1%) had the highest accuracy, while DL-G35 and DL-G50 had different degrees of overestimation (Figure 4b,e). For DRS, Poly-G35 (NSE = −0.27, Pbias = −1.49%) had the best recognition effect for SOS. POS identified by all the algorithms for DRS had large deviations, but the POS identified by the two threshold methods under the polynomial fitting (NSE = −0.58, Pbias = −1.39%) method was relatively small. EOS estimated by Poly-G35 (NSE = 0.29, Pbias = −0.61%) for DRS was in good agreement with the results of the station GPP, while the deviation of the DL-CUM results was the largest (Figure 4i,l).

In summary, among the five methods based on EVI, although DL-CUM achieved the best performance in identifying the phenological parameters of TYG (Figure 4g–i), it was not suitable for DRS. Poly-35 showed the strongest potential when based on EVI data in DRS (Figure 4j–l). In addition, among the three phenological parameters, SOS and POS had the highest consistency with real SOS, while EOS had the highest uncertainty.

For SIF data (Figure 5), there was also a larger error for DRS than for TYG. For TYG, Poly-G50 (NSE = 0.12, Pbias = 0.22%) had the highest consistency in identifying SOS, while DL-CUM and Poly-G35 had large degrees of overestimation (Figure 5g,m). DL fitting combined with 35% and 50% thresholds (NSE = −0.01, Pbias = −0.06%) gave the best POS results, followed by DL-CUM; polynomial fitting was the worst. The Poly-G35 (NSE = 0.02, Pbias = 0.08%) method has the highest accuracy in identifying EOS in the TYG region (Figure 4l). For DRS, all methods overestimated the SOS, and DL-G50 (NSE = −3.43, Pbias = 5.99%) performed relatively well (Figure 4f). In POS estimation, although all methods had different degrees of underestimation, the POS recognition results of DL fitting with the 35% and 50% thresholds were closest to the real POS (NSE = −1.0, Pbias = −2.92%). In terms of EOS, all methods showed overestimation; however, DL-CUM performed the best (NSE = −2.72, Pbias = 2.71%). In summary, all methods could well determine SOS and POS, but the EOS results showed significant uncertainty.

The optimal data sources for the three phenological parameters in TYG were all SIF, while those for DRS were all EVI. The optimal identification scheme for SOS in TYG coverage areas was Poly-G50 based on SIF data (NSE = 0.12, Pbias = 0.22%), and the optimal combination for identifying POS was the polynomial fitting method with SIF data (NSE = −0.01, Pbias = −0.06%). The optimal identification scheme for EOS was Poly-G35 based on SIF data (NSE = 0.02, Pbias = 0.08%). SOS for DRS areas had the highest accuracy using Poly-G35 with the EVI data (NSE = −0.27, Pbias = −1.49%). The optimal solution for POS was DL fitting with the EVI data (NSE = −0.58, Pbias = −1.39%). For EOS, the combination of EVI data and Poly-G35 had the highest accuracy (NSE = 0.29, Pbias = −0.21%). In summary, the potential of SIF for estimating phenology can only be realized in areas with relatively large vegetation coverage, with EVI-based phenology more suitable for DRS with sparse vegetation. The polynomial fitting method performed better than the DL fitting for both vegetation types.

3.2. Phenology Estimated from EVI

The differences among the five phenology algorithms based on EVI data were compared across the whole study area. SOS recognized by Poly-G35 (DRS) and Poly-G50 (TYG) was closest to the real SOS and was defined as SOS_true. From Figure 6d, although DL-G50 overestimated SOS in central and northwestern regions, it was more consistent in regions where SOS_true exceeded 70%. Meanwhile, in Figure 6g, SOS_TYG from the DL-CUM method differed from that of Poly-G35 by up to 28.48 days; the smallest difference was for FV (15 days). The difference between SOS_DL-G35 and SOS_true for DRS was at least 1 day (Figure 7a). Similarly, in terms of POS estimation, except for DL-CUM, the other methods were defined as POS_true; DL-CUM underestimated POS in the entire area. Figure 6h shows that POS in most regions of the MP was earlier than 185 days (Figure 7b). Under the same threshold, POS based on DL fitting in the central and northeastern MP was earlier than that based on polynomial fitting. In terms of EOS, DL-CUM (TYG) and Poly-G35 (DRS) were defined as EOS_true under the corresponding vegetation coverage types, and EOS estimated by Poly-G35 had the largest difference from EOS_true (Figure 6i,l). In particular, TYG had the largest difference. The earliest EOS detected by DL-CUM was 249 days, 51 days earlier than Poly-G35. In FV, the EOS of DL-G35 was different by 33 days from that of DL-CUM (Figure 7c). Under the same threshold, the EOS of the MP eastern FV based on DL fitting was later than that of polynomial fitting, showing the opposite characteristics of the peak period.

3.3. Phenology Estimated from SIF

The differences in the five phenology algorithms based on SIF data were compared across the whole study area. SOS recognized by DL-G50 (DRS) and Poly-G50 (TYG) was closest to the real SOS and defined as SOS_true. As shown in Figure 8g, the differences between the earliest SOS and SOS_true of DL-CUM in central and western regions of the MP were the largest; for DRS, it reached 33.41 days. As vegetation cover increased, the differences among the different methods reduced (Figure 9a). SOS_DL was later than SOS_Poly (Figure 8b,e), and the 50% threshold gave a later estimate than the 35% threshold (Figure 8a,d).

POS identified by the DL fitting method was defined as POS_true. From Figure 8b,e,h,k,n, POS_Poly was much later than POS_DL. In particular, DRS in the central and western parts of the MP had the largest differences, reaching 13 days (Figure 9b). Only the DL-CUM method produced a difference of ~10 days compared with other methods in >40% of the regions. This indicates that all phenology identification algorithms have high accuracy in estimating POS but that DL-CUM has the most limitations.

DL-CUM (DRS) and Poly-G35 (TYG) were defined as EOS_true under the corresponding vegetation coverage types. From Figure 9c, the difference was largest in FV, and Poly-G35 EOS was > 27 days later than EOS_true. From Figure 8f,i, >70% of the DL-G50 and DL-CUM areas had relatively consistent spatial distributions, which confirms that DL-G50 may be suitable for identifying phenology in such areas. All other methods differed from the selected EOStrue. In particular, Poly-G35, while achieving the best performance in validation with flux towers in TYG, had the lowest agreement with other methods in other vegetation types. This also shows that EOS always has greater uncertainty than SOS and POS in the identification scheme.

4. Discussion

4.1. Algorithm and Data Performance

Since the vegetation morphological phenology represented by EVI and the physiological phenology of vegetation photosynthesis represented by SIF correspond to different attributes of vegetation [40], there remain large differences in the identification of phenological parameters between different data sources [14,41]. In comparing the phenological identification results of EVI and SIF with those of GPP at the flux tower site, we found that the phenological identification accuracy of SIF was higher overall. This difference may be attributed to the close association between SIF and GPP. There is a certain hysteresis effect between leaf growth and photosynthesis [42]. Differences in phenological accuracy between SIF and EVI are advantageous in different vegetation types and more obvious in grasslands with sparse vegetation. Specifically, EVI is superior to SIF in terms of phenology identification in desert steppe areas (Figure 4 and Figure 5), which may be due to SIF being a weak signal and having low values and great uncertainty for sparse vegetation areas. EVI, on the other hand, has a wider application prospect in the identification of phenological parameters in arid, sparsely vegetated areas [11,43]. The value of SIF in phenology in areas with better vegetation coverage is currently being affirmed [7,44]. When estimating the phenology of typical grasslands, the performance of SIF data in typical grasslands is superior, suggesting that, under certain land cover conditions, SIF can more effectively remove the interference of ground background nonfluorescent targets compared with EVI, proving SIF has greater potential than vegetation-greenness-based indices in large-scale phenological simulation studies [45,46] and showing the limitations of SIF in estimating phenology in arid desert steppes.

Comparing the different fitting methods, the phenological identification results of the two sets of data showed that the Poly fitting method was closer to the phenological identification results of the GPP than DL fitting methods (Figure 4 and Figure 5), consistent with the results of Snyder et al. [47]. After Poly fitting, the phenological results of 35% and 50% thresholds were compared, showing that the accuracy of the 35% threshold was higher, possibly due to the time span of 50% dynamic threshold identification being short, thereby overestimating SOS and underestimating EOS. On the contrary, 50% after DL fitting was more accurate than 35%, possibly because the shape of the plant-growth curve fitted by DL was narrower, which does not match the real phenological changes, and the phenology cannot be accurately estimated, especially in the desert steppe (Figure 6 and Figure 8). At the same time, our findings demonstrated that the Poly fitting method has higher applicability in arid/semiarid areas, especially in arid desert steppe areas [48]. However, compared with the phenological identification accuracy in typical grasslands, the identification accuracy of phenology in the desert steppe was lower, due to the smaller number of years of stations available in the desert steppe (Table 1) and larger fitting errors [49].

Assuming that the difference caused by the number of stations is low, compared with the phenological parameters identified by GPP, the EOS errors generated by different phenological algorithms were usually larger than those for SOS and POS, and EOS had more uncertainty than SOS, consistent with previous findings [50,51]. It could be attributed to the fact that seasonal decoupling of physiological states and green information mainly occurs in autumn and that autumn vegetation has a more sensitive response to climate [52]. The phenology recognition accuracy of different phenology recognition algorithms is significantly different due to differences in algorithms and parameters in different fitting methods [30]. The Poly fitting method can better track the growth trajectory curve of vegetation than the DL fitting method (Figure 4 and Figure 5). The vegetation growth cycle fitted by DL is shorter than that of Poly, since the SIF curves of the same vegetation type are more “steep and sharp”, while the EVI curves are more “slow and wide”. Therefore, different thresholds have less influence on SIF phenology and greater influence on EVI phenology.

4.2. Limitations and Prospects

An in-depth analysis of the differences in the identification of phenology arising from using different data and methods was conducted in this study; however, there remain some deficiencies that must be addressed in future research. The use of a data source with a coarse resolution of 0.05° indicates that the phenological differences of mixed species in a 5 × 5 km pixel are negligible. Since the plant function types of different species may vary, a phenological value was used to represent a 25 km² area, which was not highly representative. Therefore, future studies should combine high-resolution data such as Sentinel-2 to improve the estimation accuracy of LSP [10]. In addition, EVI data provided every 16 days may miss important nodes of vegetation phenology changes, which may explain why it was different from the phenology identified using the same algorithm for SIF data with 8-day resolution [53]. To capture phenological information more accurately, data with shorter time steps should be used in the future.

Secondly, due to the limitation of the number and range of flux stations in the study area, only the 15-year data from the three stations of TYG and DRS in Inner Mongolia were used for verification, which has two shortcomings: firstly, the typical grassland area was vast, and having only two sites was insufficient to represent the typical grassland. Secondly, due to the lack of verification of large-scale FV and MEG, the applicability of data and methods cannot be evaluated, and the optimal data or algorithms found following the verification of existing stations do not correspond to their applicability to the Mongolian Plateau as a whole. To overcome these limitations, future research should collect and record more field measurement data for better reliability.

In future phenology research, it is necessary to combine multiple remote sensing data sources at finer time and space scales, make full use of the advantages of SIF and EVI in identifying phenology in different vegetation coverage areas, and comprehensively consider vegetation physiology. Due to the difference between physiological and morphological phenologies, research on phenological parameters over large spatial ranges is another future direction [54].

5. Conclusions

Using two types of remote sensing data combined with five phenological fitting methods to identify phenology, this study conducted a difference analysis of phenology identification and screened the optimal combination of phenology data sources and algorithms for the MP. The main conclusions are as follows:

(1)

The optimal methods for identifying the SOS, POS, and EOS of typical grassland areas were Poly-G50 (NSE = 0.12, Pbias = 0.22%), DL-G35/50 (NSE = −0.01, Pbias = −0.06%), and Poly-G35 (NSE = 0.02, Pbias = 0.08%), respectively, based on SIF data. The best methods for identifying the SOS, POS, and EOS of desert steppe areas were Poly-G35 (NSE = −0.27, Pbias = −1.49%), Poly-G35/50 (NSE = −0.58, Pbias = −1.39%), and Poly-G35 (NSE = 0.29, Pbias = −0.61%), respectively, based on EVI data.

(2)

The data source was the main contributor to differences in estimates of vegetation phenology on the MP. EVI data were more compatible than SIF data for identifying phenology in sparsely vegetated areas.

(3)

Satellite remote sensing can accurately capture phenological changes, but different algorithms will affect the recognition accuracy. In the same data, the sparser the vegetation, the greater the error of the phenology identification scheme.

The results of this study make an important contribution to improving current phenology estimation models, establishing a phenology estimation system suitable for the MP, and perfecting vegetation phenology observations.