Comparing Different Spatial Resolutions and Indices for Retrieving Land Surface Phenology for Deciduous Broadleaf Forests

Abstract

Deciduous broadleaf forests (DBF) are an extremely widespread vegetation type in the global ecosystem and an indicator of global environmental change; thus, they require accurate phenological monitoring. However, there is still a lack of systematic understanding of the sensitivity of phenological retrievals for DBF in terms of different spatial resolution data and proxy indices. In this study, 79 globally distributed DBF PhenoCam Network sites (total 314 site-years, 2013–2018) were used as the reference data (based on green chromaticity coordinates, GCC). Different spatial resolutions (30 m Landsat and Sentinel-2 data, and 500 m MCD43A4 data) and satellite remote sensing vegetation indices (normalized difference vegetation index, NDVI; enhanced vegetation index, EVI; and near-infrared reflectance of vegetation, NIR_V) were compared to find the most suitable data and indices for DBF phenological retrievals. The results showed that: (1) for different spatial resolutions, both 30 m Landsat–Sentinel-2 data and 500 m MODIS data accurately captured (R² > 0.8) DBF phenological metrics (i.e., the start of the growing season, SOS, and the end of the growing season, EOS), which are associated with the comparatively homogeneous landscape pattern of DBF; (2) for SOS, the NIRv index was closer to GCC than EVI and NDVI, and it showed a slight advantage over EVI and a significant advantage over NDVI. However, for EOS, NDVI performed best, outperforming EVI and NIRv; and (3) for different phenological metrics, the 30 m data showed a significant advantage for detecting SOS relative to the 500 m data, while the 500 m MCD43A4 outperformed the 30 m data for EOS. This was because of the differences between the wavebands used for GCC and for the satellite remote sensing vegetation indices calculations, as well as the different sensitivity of spatial resolution data to bare soil. This study provides a reference for preferred data and indices for broad scale accurate monitoring of DBF phenology.

1. Introduction

Vegetation phenology is a sensitive indicator of global environmental change, and the annual seasonal cycle of vegetation growth and senescence is used as an indicator for environmental change assessment [1,2]. Deciduous broadleaf forests (DBF), as a globally distributed vegetation biome, play a prominent role in the terrestrial response to global environmental change [3,4,5]. For example, vegetation phenology changes show a positive correlation with the distance to high-temperature areas (e.g., cities), indicating that vegetation phenology has a significant response to changes in the surrounding environment [6]. In addition, it has been demonstrated that global environmental change can be accurately mapped onto changes in the growth and development cycle of vegetation at a broad scale [7]. Therefore, achieving accurate monitoring of DBF at a broad scale is important for assessing environmental dynamics.

Early monitoring of phenology relied mainly on manual observation records, which were abundant in detail but difficult to observe on a broad scale; also, subjective intention and a priori knowledge could affect the phenology results to different degrees [8,9]. In recent years, with the development of remote sensing, phenological observation has entered a new stage, and data are collected by various sensors rather than by humans [10]. Methods for near-surface monitoring based on phenology cameras [11] and remote sensing monitoring based on satellite sensors, such as land surface phenology (LSP) [12,13], are increasingly applied for phenological monitoring.

Current satellite remote sensing-based LSP monitoring mainly uses 500 m MODIS data [14,15,16] and the traditional normalized difference vegetation index (NDVI) [17] and enhanced vegetation index (EVI) [18], whereas, to quantify phenological results, the start of the growing season (SOS) and end of the growing season (EOS) [19,20] phenology metrics are mostly used. There are numerous studies in this field, and the main conclusions are that LSP based on satellite remote sensing can detect the extension of the growing season, advance of SOS, and lag of EOS [21,22]. However, different studies have pointed out that, based on different data sources, the phenological results obtained by different indicators are quite variable, and can even reach completely opposite conclusions. For example, Geli Zhang et al. [23] noted in their study that SOS was not always advanced or lagged, but rather exhibited fluctuations, while the results from different data sources showed significant differences. Another study highlighted the effect of temporal resolution on SOS in DBF, with low temporal resolution data detecting false SOS advances [24]. Recent research indicated that differences in EOS may be due to different vegetation indices and algorithms [25]. The EVI generally performs better than the NDVI for these measurements and shows high accuracy; however, the emergence of the new index may challenge this status [26,27]. For example, in the near-infrared reflectance of vegetation (NIRv) index, the systematic evaluation of NIRv retrieval of vegetation phenology is not sufficient yet; however, it has been shown that NIRv can accurately estimate vegetation phenology and can estimate vegetation phenology more effectively than EVI and NDVI [28,29].

In addition, with the advent of Landsat and Sentinel-2 data, higher spatial resolution data have been applied to the phenology study of DBF. For example, Melaas E K et al. [30] conducted a study based on Landsat TM/ETM+ data to examine interannual variation in DBF phenology and the results were in excellent agreement with the field measurements. Although Landsat can provide observations at 30 m resolution, its 16-day revisit period limits the quality of the time series and is more susceptible to cloud and rain [31]. In a recent study, Sentinel-2 data were applied to monitor vegetation phenology in Nordic forests [32]. Their high temporal–spatial resolution (10 m, 5-day) data can provide a richer, more detailed portrayal to discern landscape heterogeneity [33]. However, Sentinel-2 has few historical archives and cannot perform long time series analysis. The fusion of the two data sources provides new opportunities: the fused Landsat and Sentinel-2 data can provide global observations with an average revisit period of three days, which are more conducive to the spatial and temporal analysis of DBF phenology [34]. This compensates for Landsat’s revisit period shortcomings, and the physical accuracy is greatly improved compared with a single data source [35].

Near-surface remote sensing observation provides new ideas for DBF phenology monitoring. Previous studies of vegetation phenology retrieval based on satellites have used near-surface observations as validation samples, including manual observations and repeated observations by digital cameras; however, manual observations have been gradually phased out in favor of more efficient and objective digital cameras (e.g., PhenoCam network) [8,9]. Nagai S et al. [36] used digital cameras to monitor canopy phenology in DBF, and the frequent high-quality observations of digital cameras gave remote sensing phenology results that were close to real-time changes. It has been demonstrated that near-surface remote sensing data can play a useful role in validating vegetation phenology indices obtained from satellite remote sensing data [37,38]. Several studies on vegetation phenology retrieval by satellite have been conducted using PhenoCam network observations as ground verification, and their objectivity is increasingly recognized [11,39,40,41]. However, it is difficult to support a broad scale study with PhenoCam data alone, and previous studies have applied PhenoCam data more to validating satellite phenology [37,40,42] and checking phenology models [39,43].

Although all the above studies achieved DBF phenology monitoring, the different spatial resolutions (mixed pixels) and validation samples of remote sensing images can lead to bias in vegetation phenology retrieval results. There is a lack of comprehensive and systematic assessment of whether DBF phenology at 30 m resolution is comparable with that at 500 m; whether remote sending monitoring is consistent with the PhenoCam phenology results; and which index better characterizes DBF phenology in different regions of the world. Therefore, to address the above issues, the objectives of this study were to: (1) investigate how different spatial resolution data (30 m spatial resolution Landsat and Sentinel-2 fusion data, 500 m spatial resolution MODIS data) differ in the retrieval of DBF phenology over a broad scale, and whether the results are consistent with those of near-surface remote sensing (PhenoCam); (2) compare the results of the inversion of DBF phenology based on different vegetation indices (EVI, NDVI, and, NIRv), especially whether NIRv can perform well in this regard compared with the more commonly applied NDVI and EVI; and (3) investigate the differences in DBF phenology with different phenological metrics (SOS and EOS) and their driving principles. Our study provides a new perspective for studying DBF phenology at a finer scale and helps to compare and assess the spatial and temporal patterns of vegetation phenology across scales in the context of global change.

2. Materials and Methods

2.1. Study Area

The study area was defined by the monitoring range of 81 PhenoCam network sites, which were widely distributed among woodlands between 30°N and 55°N and 125°W and 0° (Figure 1,

Table 1. Information on deciduous broadleaf forests (DBF) sites in the PhenoCam network.

Site	Latitude	Longitude	Site-Years 1
shiningrock	−82.774967	35.390159	14
smokypurchase	−83.0733	35.586	14
smokylook	−83.943113	35.632529	14
mammothcave	−86.10194444	37.18583333	13
nationalcapital	−77.069496	38.88818	13
boundarywaters	−91.495506	47.946702	13
caryinstitute	−73.7341	41.7839	11
harvard	−72.1715	42.5378	11
bartlettir	−71.2881	44.0646	11
acadia	−68.26083333	44.37694444	11
dollysods	−79.427041	39.099529	10
morganmonroe	−86.4131	39.3231	10
proctor	−72.866	44.525	10
queens	−76.324	44.565	10
umichbiological2	−84.6976	45.5625	10
joycekilmer	−83.795	35.257	9
howland2	−68.7418	45.2128	9
umichbiological	−84.71382	45.55984	9
woodshole	−70.6432	41.5495	8
hubbardbrook	−71.701	43.9438	8
snakerivermn	−93.24467	46.120556	8
drippingsprings	−116.8	33.3	7
oakridge2	−84.3323	35.9311	7
tonzi	−120.9658861	38.43091667	7
missouriozarks	−92.2	38.7441	7
northattleboroma	−71.31056	41.98369	7
bostoncommon	−71.064145	42.355912	7
readingma	−71.1272	42.5304	7
harvardlph	−72.185	42.542	7
turkeypointdbf	−80.5576	42.6353	7
willowcreek	−90.07912	45.805986	7
coweeta	−83.4275	35.0592	6
alligatorriver	−75.9038	35.7879	6
springfieldma	−72.585972	42.135162	6
ashburnham	−71.926	42.6029	6
downerwoods	−87.88076	43.07938	6
uwmfieldsta	−88.0229	43.38709	6
hubbardbrooksfws	−71.7407	43.9269	6
russellsage	−91.974322	32.456961	5
dukehw	−79.100371	35.973583	5
bullshoals	−93.06663	36.562833	5
shalehillsczo	−77.9041	40.6658	5
worcester	−71.8428	42.2697	5
arbutuslake	−74.23322	43.98207	5
laurentides	−74.0055	45.9881	5
canadaOA	−106.19779	53.62889	5
silaslittle	−74.596	39.9137	4
bbc2	−72.185	42.542	4
sanford	−84.464482	42.726777	4
lacclair	−71.66957778	46.95209167	4
marcell	−93.46925	47.5139	4
NEON.D02.SERC.DP1.00033	−76.56001	38.89008	3
bbc5	−70.6432	41.5495	3
bbc1	−72.174359	42.53508	3
tfforest	−70.9505	43.1086	3
arbutuslakeinlet	−74.24527	43.99336	3
bbc7	−71.2881	44.0646	3
NEON.D08.LENO.DP1.00033	−88.16122	31.85388	2
NEON.D08.DELA.DP1.00033	−87.803877	32.541727	2
NEON.D11.CLBJ.DP1.00033	−97.57	33.40123	2
NEON.D07.GRSM.DP1.00033	−83.50195	35.68896	2
NEON.D07.ORNL.DP1.00033	−84.282588	35.964128	2
asuhighlands	−81.7032	36.2076	2
NEON.D02.SCBI.DP1.00033	−78.139494	38.892925	2
usgsreston	−77.3676	38.9471	2
NEON.D02.BLAN.DP1.00033	−78.041788	39.033698	2
macleish	−72.6804	42.4484	2
NEON.D01.HARV.DP1.00033	−72.17265	42.536911	2
NEON.D01.BART.DP1.00033	−71.287375	44.063869	2
willamettepoplar	−123.1823	44.1368	2
NEON.D05.TREE.DP1.00033	−89.58572	45.49373	2
NEON.D05.UNDE.DP1.00033	−89.537254	46.23391	2
millhaft	−2.29883	52.800796	2
donanapajarera	−6.4432	36.9962	1
NEON.D07.MLBS.DP1.00033	−80.524847	37.378314	1
robinson2	−83.1576	37.4671	1
pace	−78.2739	37.9229	1
columbiamissouri	−92.1997	38.7441	1
greenridge1	−78.4067	39.6905	1
hubbardbrooknfws	−71.7762	42.958	1
NEON.D05.STEI.DP1.00033	−89.58637	45.50894	1

¹ Site-years represents the number of years of data from 2003 to 2018 used in this study.

). All but two sites (“millhaft” and “donanapajarera” in Western Europe) were established on the southeast coast, east-central, and the southwest coast of North America. Most of the sites were in a temperate continental climate, with cold winters and warm summers, and the forests were homogeneous DBF with obvious seasonal changes in phenology [44]. The southeastern sites were in a subtropical humid climate, mostly influenced by monsoons, and a small amount of evergreen broadleaf forests could be mixed in the forests [45]. The northern sites could be mixed with a small number of coniferous species.

2.2. PhenoCam Data

In this study, the DBF phenological periods from near-surface remote sensing observations captured by digital cameras were used as a validation for the satellite remote sensing phenological retrieval results. Specifically, we chose the relevant data from the PhenoCam Data (v2.0) (https://PhenoCam.nau.edu/PhenoCam_explorer/; last accessed on 7 September 2022) [39]. The data are collected from high-resolution cameras set up at over 700 sites in the PhenoCam network. The sites are typically elevated and have a fixed field of view (FOV) for continuous observation, acquiring images with very short time intervals (up to 30 min) in the time series [46].

We selected 81 sites from which DBF could be well observed. The data acquired from these sites spanned more than 18 years (2001–2018). However, considering the lower observation quality and a smaller quantity of early PhenoCam data and satellite data, the data from 2013 to 2018 were mainly used for a total of 79 sites. Two sites (“drippingsprings” and “columbiamissouri”) were dropped due to the lack of data from 2013 to 2018.

In the processing of the PhenoCam data, the PhenoCam team officially set the region of interest (ROI) for each image to exclude the influence of non-canopy imaging elements on the results. The ROI is used to extract vegetation indices, such as the green chromatic coordinate (GCC) index, and then construct 1-day or 3-day GCC time series. The team also calculated the phenology periods such as greenness rise (SOS) and fall (EOS) of greenness transition with different amplitude thresholds based on the time series. It is worth noting that the location information of the PhenoCam phenology cameras provided on the official website did not match the location information of the actual observation area of the camera (the official region of interest (ROI) provided). Therefore, we also corrected the location information of each station according to the observation direction of the phenology camera and referred to the images taken by the PhenoCam camera to provide satellite image calibration using location information closer to the observation area of the PhenoCam camera. In contrast to the satellite remote sensing images, the images acquired by PhenoCam were mainly sensitive to the visible band; thus, the vegetation index used for the PhenoCam data was based on the RGB three-color channel calculation in the visible band. The $GCC$ index can provide a good indication of the greenness of the canopy layer in the ROI by calculating the proportion of green digital numbers to the total red, green, and blue digital numbers in the ROI, which, in turn, can reflect the change in DBF phenology [47]. The data derived from PhenoCam in this study are the phenology data retrieved at 10%, 25%, and 50% thresholds based on 3-day time series compounded by 90th percentile $GCC$ values, and the $GCC$ calculation formula is as follows:

$$GCC=\frac{{Green}_{DN}}{{Red}_{DN}+{Green}_{DN}+{Blue}_{DN}}$$

(1)

where ${Red}_{DN}$, ${Green}_{DN},$ and ${Blue}_{DN}$ represent the three-color RGB channels extracted from the image, and the upper and lower fractions are the digital numbers of the corresponding color channels. We can also calculate the chromaticity of other color channels using a similar formula [48].

2.3. Satellite Remote Sensing Data

2.3.1. Landsat and Sentinel-2 Data

Landsat data and Sentinel-2 data at 30 m spatial resolution were acquired from the Google Earth Engine (GEE) platform. Although the revisit period of Landsat is 16 days, the actual temporal resolution of the available data is often more than 16 days owing to the quality of the image. The Sentinel-2 thermal infrared sensor data have similar bands to Landsat8 OLI data and a shorter revisit period; thus, the fusion of the two data sets can solve the problem of long revisit times [49]. Specifically, the algorithm was built and based on GEE platform implementation. First, the Sentinel-2 images were resampled to maintain the same spatial resolution (30 × 30 m) as the Landsat image using the bicubic resampling method. Then, the two data points with the closest missing points in the time series were selected, and a complete time series was reconstructed using ordinary least squares regression coefficients (

Table 2. The transformation coefficients for Landsat and Sentinel-2 data bands.

Bands	Landsat-7	Bands	Sentinel-2
Blue	Landsat8(480 nm) 1 = 0.0003 + 0.8474 × Landsat7(485 nm) 1	Blue	Landsat8(480 nm) 1 = 0.0003 + 0.9570 × Sentinel2(490 nm) 1
Red	Landsat8(655 nm) 1 = 0.0061 + 0.9047 × Landsat7(660 nm) 1	Red	Landsat8(655 nm) 1 = 0.0041 + 0.9533 × Sentinel2(665 nm) 1
NIR	Landsat8(865 nm) 1 = 0.0412 + 0.8462 × Landsat7(835 nm) 1	NIR (Band 8A)	Landsat8(865 nm) 1 = 0.0077 + 0.9644 × Sentinel2(865 nm) 1
SWIR	Landsat8 (1610 nm) 1 = 0.0254 + 0.8937 × Landsat7(1650 nm) 1	SWIR	Landsat8(1610 nm) 1 = 0.0034 + 0.9522 × Sentinel2(1610 nm) 1

¹ The wavelengths in parentheses after Landsat7, Landsat8, and Sentinel-2 in the table are the central wavelengths of their corresponding bands.

) to transform bands 4 (red), 8A (near infrared, NIR), 11 (short-wave IR, SWIR), and 2 (blue) in the Sentinel-2 MSI data [50]. The data were divided into three-time points according to different data types: 2001–2012, with only Landsat data (excluding Landsat8), 2013–2015, with all Landsat series data, and 2015–2018, with both full series Landsat data and Landsat–Sentinel-2 fusion data.

Three vegetation indices—$NDVI$, $EVI$, and $NIRv$—were selected from the satellite remote sensing data using the following formulae:

$$NDVI=\frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+{\rho }_{RED}}$$

(2)

$$EVI=2.5\times \frac{{\rho }_{NIR}-{\rho }_{RED}}{{\rho }_{NIR}+6\times {\rho }_{RED}-7.5\times {\rho }_{BULE}+1}$$

(3)

$$NIRv=NDVI\times {\rho }_{NIR}$$

(4)

where ${\rho }_{NIR}$, ${\rho }_{RED}$, and ${\rho }_{BULE}$ represent the reflectance of the NIR band, red band, and blue band in satellite images, respectively. Both $NDVI$ and $EVI$ have been widely used for vegetation green change monitoring [51]. However, the effectiveness of the $NIRv$ index for monitoring vegetation phenology remains to be verified [28]. To ensure the spatial consistency of the extracted phenology results from satellite images with the validation samples (PhenoCam data), we set gradient buffers (30–90 m) in the GEE platform for each image according to the corrected PhenoCam site location information. Specifically, gradient buffers of 30 m, 60 m, and 90 m in radius were set as the center of the circle with the PhenoCam site location (Figure 2). Using these buffers, the NDVI, EVI, and NIRv were calculated separately and were compared with the GCC indices obtained by PhenoCam.

2.3.2. MODIS Data

The 500 m resolution MODIS data were taken from the NASA LP DAAC (LP DAAC—Data (usgs.gov); last accessed on 7 September 2022) surface reflectance product (MCD43A4 Version 6 Nadir Bidirectional Reflectance Distribution Function Adjusted Reflectance (NBAR)). The data were selected for the years spanning 2001–2018, but, to be consistent with PhenoCam data years, we determined the time range of remote sensing data to be 2013–2018. MCD43A4 are daily data, which collect the most representative data in every 16 days, with the 9th day of 16 days as the recording date. The product also provides data from vertical shots of the land surface, which can characterize canopy greenness changes and be used to construct a time series of vegetation indices. However, the MCD43A4 does not guarantee that all image elements are reliable; thus, it was necessary to use the MCD43A2 Version 6 Bidirectional Reflectance Distribution Function and Albedo (BRDF/Albedo) Quality product to perform quality assurance (QA), remove cloud, etc. for the MCD43A4 data. MCD43A2 contains the quality information for the corresponding MCD43A3 Albedo and MCD43A4 NBAR products. We filtered out the unneeded image elements based on the QA information provided by MCD43A2 for MCD43A4. Specifically, in this study, MCD43A2 QA was used, and pixels representing snow, as well as pixels with poor data quality in each band (1–5), were removed according to Quality Band pixel_qa, such as pixel_qa1, pixel_qa2, and pixel_qa3. For the processed MCD43A4 data, we then calculated the NDVI, EVI, and NIRv indices separately in Equations (2)–(4).

2.4. Phenological Retrieval for Different Data

In this study, phenological retrieval was based on the time series constructed from the vegetation indices calculated from the above data. To ensure the accuracy of phenological retrieval, we excluded the effect of missing values and outliers in the original calculation results of vegetation indices. We arranged the initial calculations of the four vegetation indices for each site–year in a time series (day of the year, DOY), filled data gaps, and performed smoothing operations on the raw data. Specifically, we first supplemented the missing values with a cubic spline interpolation and then performed the first smoothing using the Savitzky–Golay filter to initially filter out the noise caused by factors such as bad weather. A double logistic function, Equation (5), was then used to fit a double logistic (time series) curve based on the initially filtered time series and to complete the quadratic smoothing [52].

$$f\left(t\right)={\alpha }_1+\frac{{\alpha }_2}{1+e^{-{\partial }_1(t-{\beta }_1)}}+\frac{{\alpha }_3}{1+e^{-{\partial }_2\left(t-{\beta }_2\right)}}$$

(5)

where $f\left(t\right)$ represents the calculated value of the vegetation indices on DOY (t), ${\alpha }_1$ is the initial value of the function fit, which is the retrieval value in the background of the vegetation indices for DBF (generally in winter), and ${\alpha }_2$ and ${\alpha }_3$ correspond to the retrieval values of the vegetation indices for SOS and EOS, respectively. ${\partial }_1$ represents the inflection point where the curve begins to rise (SOS), and ${\partial }_2$ represents the inflection point where the curve begins to fall (EOS); ${\beta }_1$ and ${\beta }_2$ are the slope coefficients of SOS and EOS, respectively [53].

The relative threshold method was used to extract the phenological metrics (SOS and EOS) [54]. For double logistic curves, we chose 10–25–50% amplitudes to extract. The first and last DOY that exceeded our set thresholds were defined as SOS and EOS, respectively.

2.5. Comparisons of Phenological Results

For time series variation, the data from each vegetation index calculation can be lengthy and difficult to compare. For interannual variation, SOS and EOS based on time series extraction were analyzed in this study; for intra-annual variation, we used visual comparison (plotting double logistic curves for each year). For different spatial resolutions, vegetation indices, and phenological metrics, we introduced the statistical parameters root mean square error ($RMSE$) and R-square ($R^2$). The $RMSE$ and $R^2$ between Landsat–Sentinel-2/MODIS phenology results and PhenoCam phenology results were calculated Equations (6) and (7) to quantitatively compare the phenology results of the two data sets.

$$RMSE=\sqrt[2]{\frac{1}{m}{\sum }_{i=1}^m{(y_i-{\widehat{y}}_i)}^2}$$

(6)

$$R^2=1-\frac{{\sum }_i{({\widehat{y}}_i-y_i)}^2}{{\sum }_i{({\overline{y}}_i-y_i)}^2}$$

(7)

where $m$ represents the total number of stations, $i$ represents the ith station, and $(y_i-{\widehat{y}}_i)$ represents the difference between the phenological period retrieved from the satellite remote sensing and the phenological period retrieved from PhenoCam. $RMSE$ takes the value [0, +∞), where smaller values showed that the satellite phenology results were closer to PhenoCam and had less dispersion; $R^2$ takes the value [0, 1], whereby the closer the fit result was to 1, the better the satellite phenology results matched PhenoCam.

3. Results

3.1. Time Series Variation of Satellite Vegetation Indices

The time series curves of each site were prepared to extract the phenology. As shown in Figure 3, we fitted the time series curve in 6 typical sites using the double logistic curve (2018) to extract SOS and EOS at 30 m and 500 m spatial resolution. The results showed that the 30 m spatial resolution data was more lagged in SOS and EOS relative to the 500 m spatial resolution data in each index. Combined with Figure 4, our extracted phenological period (6 typical sites, 2018) broadly matched the temporal changes in the site photographs captured by PhenoCam.

In Figure 5, a mixed time series of all sites from 2013 to 2018 is shown. The inversion values of the vegetation indices exhibited a larger amplitude in the Y-axis direction and a more uniform distribution. The spatial variability brought by the latitude of each site (Figure 6) was weakened, and the time series variation of the overall study area was highlighted. The mixed time series shows more clearly that the time series alignment of each index for the 30 m data was coarser than that of the 500 m data (MODIS), and the MODIS time series was smoother. The remote sensing indices (GCC, NDVI, EVI, NIRv) had similar characteristics to the time series, showing a slow increase (February–March) followed by a fast increase (April–June) and then a fast decrease (July–October) followed by a slow decrease (October–November) within the year. The interannual variation was insignificant for each index, indicating that our results were stable, and the method was reliable. The mixed site time series clearly showed that SOS and EOS were more lagged at 30 m compared with 500 m, which suggested that all of our site year phenology results supported the conclusion that 30 m was closer to PhenoCam relative to 500 m for SOS.

3.2. Phenological Comparisons at Different Spatial Resolutions

We processed data from all remote sensing data for 314 site-years and extracted the results of the start of SOS and EOS from time series fitted curves in 3 gradient buffers (30–90 m) according to the optimal threshold (25%) indicated by our data results. From the overall results, the phenology results of different indices in different data showed quite obvious spatial variability, with SOS differences of up to 92–114 days and EOS differences of up to 74–149 days between different stations. The SOS increased nonlinearly with increasing latitude, and the EOS decreased nonlinearly with increasing latitude, but there was no obvious linear change with longitude. We speculated that the changes were mainly influenced by the difference in solar radiation between different latitudes (Figure 6).

To visually demonstrate the difference in sensitivity between the 30 m and 500 m spatial resolution data for the detection of phenology in broad scale DBF, we performed a linear fit of the phenological periods extracted from the GCC and each satellite remote sensing index. The results indicated that both 30 m and 500 m scales performed well (as shown in Figure 7). Among the 36 control groups [satellite (30 m, 500 m) phenology results and PhenoCam phenology results] in the full gradient buffer, 26 groups reached R² above 0.8, which was related to the relatively homogeneous landscape pattern of DBF. The highest value of R² = 0.87 occurred in the buffer_90 m SOS-30 m spatial resolution-NIRv and GCC control (Figure 7), and the lowest value of R² = 0.69 occurred in the buffer_60 m EOS-30 m spatial resolution-NIRv and GCC control (Figure 7).

On the whole, the phenology results at 30 m differed substantially from those at 500 m, and it is worth noting that the SOS results were significantly better than the EOS results. For the SOS, the 30 m effect on each remote sensing vegetation index (NDVI, EVI, NIRv) was generally better than at 500 m, and was not affected by the selection of the vegetation index and buffer gradient. All groups of 30 m (the average R² was above 0.84) fit better than 500 m (the average R² was above 0.81) in the nine control groups (30 m and 500 m) containing all the indices. In summary (

Table 3. Accuracy comparison of phenology results with different spatial resolutions (30 m and 500 m) based on PhenoCam GCC.

Buffer	Metric	PhenoCam_GCC	Landsat and Sentinel-2			MODIS
			EVI	NDVI	NIRv	EVI	NDVI	NIRv
30 m	SOS	1	0.87	0.79	0.87	0.84	0.75	0.85
	EOS	1	0.77	0.85	0.71	0.85	0.85	0.85
60 m	SOS	1	0.86	0.79	0.87	0.84	0.75	0.85
	EOS	1	0.80	0.85	0.69	0.85	0.85	0.85
90 m	SOS	1	0.86	0.79	0.87	0.84	0.75	0.85
	EOS	1	0.80	0.82	0.74	0.85	0.85	0.85

), the 500 m results lost an average of 3% of accuracy in the 30 m–90 m buffer compared with the 30 m indices. Specifically, for each site-year data, the maximum difference between the 30 m and 500 m phenological period in the same site-year could be 28 days (buffer30 m_site:bostoncommon_2018_EVI), which compared with the GCC index (DOY = 122) in the same site-year. The SOS (DOY = 121) of 30 m was closer to the GCC value, and the index of 30 m generally had a smaller and earlier difference. We also noted that both 30 m and 500 m phenological results were significantly earlier relative to the phenological period obtained by the PhenoCam phenological camera (Figure 7). This may be because the GCC index used the visible band, which was not sensitive to early bud (not yet transformed to green, Figure 8a growth.

For EOS, the 500 m result was rather better than that of the 30 m data. The fit accuracy (the average R² was above 0.85) for the 500 m was greater than or equal to 30 m (the average R² was only above 0.78) in all nine control groups (30 m and 500 m). As shown in Table 3, the 30 m phenology results lost approximately 7% accuracy on average over the 30 m–90 m buffer compared with the 500 m results. The comparison results were not changed by the buffer gradient, but the 30 m performance in EOS was worse than the 500 m performance. In addition, the phenology results based on NDVI were significantly better. In contrast to the SOS results, the satellite remote sensing-based EOS results were significantly later than the PhenoCam results (Figure 7). This demonstrated that our speculation regarding the lagging SOS results for the GCC index was convincing because the imaging in the visible band failed to identify leaves that were browned but still viable Figure 8b,c.

3.3. Differences in GCC and Satellite Vegetation Indices for Different Phenological Metrics

Figure 9 shows that the comparison of NDVI, EVI, NIRv, and GCC highlighted significant differences in phenological metrics (SOS and EOS). The performance of EVI (average R² = 0.85) and NIRv (average R² = 0.86) for SOS was significantly better than that of NDVI (average R² = 0.77) (i.e., closer to the GCC index), and NIRv showed a slight advantage over EVI, with an average accuracy improvement of 1% in the 30 m–90 m buffer. The accuracy of the NDVI-based EOS fit (average R² = 0.84) was improved significantly, better than EVI (average R² = 0.82) and NIRv (average R² = 0.78), and EVI performed better than NIRv. We also noted that the vegetation indices performed better in the SOS metric than the EOS metric, which was consistent with the view of Yang et al. [55] that EVI and NIRv were more suitable for SOS inversion. The SOS metric was better suited at 30 m spatial resolution, while the EOS did not highlight the advantages of 30 m data well.

4. Discussion

4.1. Applicability of Different Resolution Data (30 m and 500 m) for SOS and EOS

DBF are widely distributed across the Earth, and their phenological changes are good indicators of the changing global ecological environment. To compare the sensitivity of 30 m and 500 m spatial resolution data for measuring broad scale phenological changes in DBF, this study used Landsat–Sentinel-2 (30 m) and MODIS (500 m) data to retrieve the phenological periods at 79 DBF sites in the PhenoCam network. The accuracy of the results was examined using PhenoCam GCC-retrieved phenological periods as validation samples. The results showed that a good performance was achieved for both spatial resolutions, but there was a large difference in the performance of SOS and EOS.

Landsat–Sentinel-2 data phenological periods retrieved for the SOS in this study were better than the results of MODIS data, which highlighted the advantages of high spatial resolution data in DBF phenological retrieval. Increasing sensor spatial resolution can separate vegetation from other surface objects when calculating vegetation indices, and many studies have used high-resolution images to monitor vegetation phenology in previous studies [30,35,56]. Melaas E K et al. (2013) indicated that spatial resolution had a non-negligible impact on DBF phenology monitoring [30]; Vrieling A et al. (2017) suggested that the trend in phenology monitoring is to use higher spatial resolution (30 m) data [57].

However, our results also showed that for EOS, MODIS data retrieved phenological periods with an improvement in accuracy of 7% relative to Landsat–Sentinel-2 data. Many previous studies have shown that increased spatial resolution helps to retrieve more accurate phenological periods [58,59]. Our results contradict this, which may be mainly the result of the action of two factors: (1) the PhenoCam GCC index was calculated based on the visible band while the NIR band was used in all of the satellite remote sensing vegetation index calculations [48]. The involvement of the NIR band could invert the physiological changes (e.g., chlorophyll content) of the vegetation leaves, while the GCC index tended to reflect the changes in leaf greenness [25]. This may result in the GCC retrieval value of the EOS being more advanced. In fact, the growing season of vegetation obtained from satellite remote sensing observations was longer when compared with that obtained from camera observations [60,61]. (2) The spatial heterogeneity of DBF is often overlooked, and, for EOS, bare soil reflected in remote sensing images would show colors similar to browning leaves (Figure 10). The 500 m spatial resolution data covered a large area per image element, and the bare soil part would be mixed with leaves to show more browning in the index calculation [62]. The 30 m spatial resolution was significantly better than MODIS data for detail portrayal, and the calculation of indices could distinguish the image elements of bare soil mixed with leaves in the image, thus showing a lesser amount of browning [63]. For these reasons, the 500 m phenology results in the EOS might be closer to the GCC phenology results (better fit) than the 30 m results, while the 30 m phenology showed a lag compared with both GCC and 500 m phenology (slightly worse fit than the MODIS data).

4.2. Gradient Buffers: Better Phenological Results?

This study verified the results of satellite remote sensing retrieval of DBF phenological changes based on near-surface phenological remote sensing (PhenoCam site observation data). We found that the 30 m results were not substantially superior to the 500 m results when only using the 30 m radius buffer (Figure 7, EOS 9 group control). In previous studies, setting of gradient buffers can better respond to spatial heterogeneity [64,65]. We set a larger radius buffer to verify whether the results obtained at the 30 m radius buffer were accurate. Because a circular buffer with gradient radius for Landsat-Sentinel-2 30 m spatial resolution could encompass a different number of imaged elements in the gradient variation compared with the hybrid imaged elements of MODIS 500 m spatial resolution for 30 m spatial resolution data, a 60 m radius buffer can have three times more imaged elements compared with a 30 m radius (Figure 11b), which can reflect in the phenological variation.

In this study (Figure 7, Figure 9 and Figure 11), the 30 m spatial resolution phenology results showed weak fluctuations as the buffer radius increased, which did not affect the phenological comparison results of the two spatial resolutions. Only the SOS of NDVI was positively correlated with the radius from 30 m to 90 m in the buffer gradient variation, while the rest of the indices showed a weak increase with some fluctuations. We believe that this was still related to the relatively homogeneous landscape pattern of DBF, for which a smaller buffer zone could also yield realistic phenological results, highlighted by the advantage of 30 m spatial resolution data (Figure 11a). On the other hand, we believe that these weak fluctuations in phenology results may be due to vegetation mixing (e.g., shrubs or evergreen trees) (Figure 11c), and the weak increase in the rest of the indices laterally confirms that some “foreign bodies” (e.g., shrubs or evergreen trees) are being added to the DBF buffer as the buffer changes. We considered that setting gradient buffers in this study was a necessary step, and the results of different buffers confirmed that improved spatial resolution was helpful for DBF phenology monitoring. This result is a reference for setting buffer zones for satellite monitoring of DBF phenology.

4.3. Difficulties of Satellite Remote Sensing Monitoring of Phenology

On the one hand, satellite remote sensing for broad scale phenology monitoring often contains hybrid image elements, especially the 500 m data, due to the limitation of the spatial resolution and the heterogeneity of landscape patterns. However, the landscape heterogeneity of DBF is often overlooked, and its heterogeneity is mainly manifested in changes in the background of DBF during autumn brownout, such as soil bareness and seasonal changes in grass growth (Figure 10). On the other hand, satellite remote sensing for broad scale monitoring of phenology usually faces the challenge of how to check the results, and reliable ground samples are often difficult to obtain. This study aimed to examine the comparison of DBF phenology results at different resolutions and indices on a global scale. However, the validation sample was only available for the PhenoCam network of phenology cameras distributed in North America and two sites in Western Europe. Furthermore, at present, most remote sensing phenology monitoring studies use phenology cameras or images acquired by UAV as validation samples [66,67,68]. The difference in spatial resolution with the phenology camera data and the difference in spatial resolution among the validation sample data are often overlooked when comparing the spatial resolution of different satellite data. In the case of PhenoCam, for example, the spatial resolution of the site-years phenology camera data is not consistent, and bias in the validation sample may affect accuracy in precision assessment.

In addition, although we corrected the site location information based on the FOV directions of the PhenoCam camera, the ROI of the PhenoCam camera images did not correspond to the satellite remote sensing images. Specifically, the ROI of PhenoCam camera images is an irregular polygon that changes according to the distribution of DBF in the image, and the ROI of satellite remote sensing images is a set circular buffer. The phenological period of the non-overlapping part of the two ROIs is part of the error of the satellite remote sensing DBF phenological period accuracy assessment.

5. Conclusions

In the context of global warming today, long-time series retrieval of vegetation phenology can effectively monitor the stability of the environment and ecosystem locally or globally, and it be an important reference for preventing climate extremes [69]. DBF is broadly distributed and has obvious seasonal changes, making them excellent monitoring objects, and this study can provide a reliable methodological support. Our study compared the phenology results of 314 site-years extracted from 30 m spatial resolution (Landsat and Sentinel-2) and 500 m spatial resolution (MODIS) data at different gradients of buffers (30–90 m) based on NDVI, EVI, and NIRv indices. The PhenoCam GCC was used as a reference to compare the phenological results of 30 m and 500 m data, and the performance of different vegetation indices for SOS and EOS was also compared.

The study showed that satellite remote sensing could accurately capture the phenological changes in DBF. The 30 m resolution data were more sensitive to the phenological changes in DBF at SOS; while, at EOS, the phenological periods retrieved from 500 m resolution data were closer to the PhenoCam phenological results. Comparing NDVI, EVI, and NIRv, we found that NDVI performed relatively poorly, while EVI and NIRv both obtained high accuracy for the monitoring of DBF’s phenology. However, the results obtained by EVI were slightly worse compared with those of NIRv. There are few studies on the application of NIRv to vegetation phenology monitoring, and our comparison adds to existing knowledge. Our proposed method can provide a more accurate direction for broad scale DBF phenology monitoring.