Mapping and evaluation of NDVI trends from synthetic time series obtained by blending landsat and MODIS data around a Coalfield on the loess plateau

: The increasingly intensive and extensive coal mining activities on the Loess Plateau pose a threat to the fragile local ecosystems. Quantifying the effects of coal mining activities on environmental conditions is of great interest for restoring and managing the local ecosystems and resources. This paper generates dense NDVI (Normalized Difference Vegetation Index) time series between 2000 and 2011 at a spatial resolution of 30 generated better results (0.70 < R 2 < 0.76) than scheme 1 (0.56 < R 2 < 0.70) in this study area. Trend analysis was then performed with the synthetic dense NDVI time series and the annual maximum NDVI (NDVI max ) time series. The accuracy of these trends was evaluated by comparing to those from the corresponding MODIS time series, and it was concluded that both the trends from synthetic/MODIS NDVI dense time series and synthetic/MODIS NDVI max time series (2000–2011) were highly consistent. Compared to trends from MODIS time series, trends from synthetic time series are better able to capture fine scale vegetation changes. STARFM-generated synthetic NDVI time series could be used to quantify the effects of mining activities on vegetation, but the test areas should be selected with caution, as the trends derived from synthetic and MODIS time series may be significantly different in some areas.


Introduction
With an area of 648,700 km 2 , the Chinese Loess Plateau is home to about 108 million people [1]. It is one of the largest regions with severe soil erosion in the world [2]. The largest quantity of coal resources in China are produced in the Loess Plateau (Figure 1), and the total output is increasing. Coal mining activities (coal exploiting) accelerate soil erosion, due to the alteration of underground geological structures and land surface subsidence [3,4]. It is important to assess the effects of coal mining activities on the local environment. One solution to this issue is to compare the changes in vegetation productivity, which is recognized as one of the key indicators of environmental conditions, between mined and non-mined areas [5]. Therefore, long-term earth observation data with high spatial and high temporal resolution are of great interest for the assessment of environmental conditions in coalfields on the Loess Plateau.
Landsat TM/ETM+ (Thematic Mapper/Enhanced Thematic Mapper Plus) data with a 30-m spatial resolution seems to be the best data to capture vegetation complexity at a local scale compared to 4-m IKONOS (Space Imaging) data and 250-m MODIS (Moderate Resolution Imaging Spectroradiometer) data [6] and have been widely used to map vegetation changes over local and regional scales [7][8][9][10][11][12]. Landsat scenes from the same month or season of each year are used to assemble time series data for vegetation trend analysis [8,9]. However, this kind of time series often suffers from missing data in a certain year and false trend results that are caused by variation in observation time, due to the 16-day Landsat revisit cycle and frequent cloud contamination. In addition, Landsat data alone cannot be used for phenology analysis, because of the sparse and unbalanced distribution of the acquisition date [8], limiting its application in monitoring of long-term vegetation phenology change.
In contrast, satellite sensors, such as the Advanced Very High Resolution Radiometer (AVHRR), Système Pour l'Observation de la Terre (SPOT)-VEGETATION, MODIS and ENVISAT-Medium Resolution Imaging Spectrometer (MERIS), are able to acquire observation data more frequently. Cloud contamination can be reduced or even removed by compositing images within a certain number of continuous days, thus allowing robust trend and frequency analysis to obtain phenological information [8]. With these more frequent time series data, many vegetation trend analyses have been carried out at regional to global scales [13][14][15][16][17][18][19][20]. However, the spatial resolution (250-8,000 m) of these remote sensing data is too coarse for vegetation assessment at local or smaller scales [8,21].
With a spatial resolution of 250 m or 500 m, the MODIS images can be used for monitoring the overall changes of vegetation in a whole mining area. However, the width of a single coal face is usually 100-300 m and the length, 500-2,000 m. At such a scale, Landsat images with a spatial resolution of 30 m are suitable to capture vegetation changes, while the MODIS pixel size is too coarse. Generally, the coal cutting speed is 1-2 m per day (30-60 m per month); consequently, vegetation change monitoring at a temporal resolution better than one month is necessary. Therefore, it is of great interest to combine the spatial resolution of Landsat images and the temporal resolution of MODIS images. In addition, improvements in resolution are also highlighted in mapping land surface temperature [22][23][24].
Several data fusion methods have been proposed to gap-fill the trade-off of remote sensing sensors between spatial and temporal resolution [21,[25][26][27][28], among which the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM), proposed by Gao et al. [21], has been recognized as an efficient one [29]. The feasibility of using STARFM to blend Landsat and MODIS surface reflectance data to generate synthetic reflectance data with high spatial and high temporal resolution has been shown by several studies in various regions with different vegetation types [30][31][32][33][34][35]. Several studies modified the STARFM algorithm to better capture spatial and temporal land cover changes [28,[36][37][38]. Singh [31] generated a long-term Normalized Difference Vegetation Index (NDVI) time series of cropland in India and found that the synthetic time series adequately captured the phenological patterns over a growing season in eight separate years. Bhandari et al. [33] and Schmidt et al. [35] evaluated the ability of STARFM to generate long-term synthetic time series data to extract phenological parameters using TIMESAT (a software for analyzing time-series of satellite sensor data) [39] and obtained reasonably accurate results for different vegetation communities in Australia. However, few applications of long-term vegetation trend analysis have been carried out using the STARFM-generated synthetic NDVI time series data.
This study is expected to generate synthetic NDVI time series at a spatial resolution of 30 m using STARFM and to evaluate its capability for monitoring vegetation changes around a typical coalfield on the Loess Plateau. Landsat and MODIS NDVI images were used in a two-fold way: (1) directly as inputs for STARFM to generate synthetic 30-m resolution NDVI time series and (2) for evaluation of the accuracy of the synthetic NDVI time series produced. Synthetic 30-m spatial resolution NDVI images were evaluated against four independent Landsat scenes covering different growth stages in 2000. Time series data from 2000 to 2011 were generated based on a limited number of Landsat scenes, and trend analyses were conducted for both the synthetic dense NDVI time series and annual maximum NDVI (NDVI max ) time series. The accuracy of these high spatial resolution trends was evaluated by comparing with the corresponding trends derived from MODIS NDVI time series. The methods presented provide a framework for the analysis of trends from synthesized data using the STARFM algorithm.

Study Area
The Datong Coalfield is one of the largest coalfields in China, which is located in Shanxi Province in the eastern part of the Loess Plateau ( Figure 1). The geomorphic landscape in the western part of this region is characterized by ravines and gullies, which are typical on the Loess Plateau. In contrast, the landscape in the eastern part, where Datong city is located, is much flatter. The main vegetation types are irrigated and non-irrigated croplands, shrubs, broadleaf forest and grassland. The land cover map ( Figure 2) was created using the Support Vector Machine (SVM) method by Zhao et al. [40], having an overall accuracy of over 85%. According to this classification map, the accuracy of STARFM-generated synthetic data will be assessed among different vegetation types. For this purpose, the 85% classification accuracy is adequate. The growing season of the natural vegetation types lasts from April to November, while farming periods are much shorter, from June to August. The annual precipitation in this region is about 250-450 mm.

Figure 2.
Land cover classification map of the study area (modified/adopted from Zhao et al. [40]).
Coal mining activities in the Datong Coalfield has a long history, dating from ancient China. Mining occurs in the Jurassic and Carboniferous-Permian strata, with areas of 772 km 2 and 1,739 km 2 , respectively. Until recently, all the mining activities were located in the Jurassic strata, as shown by the boundary in Figure 1. The majority of coal resources in these strata were found about 300 m below the surface, and about 690 km 2 of land has experienced surface subsidence, due to the mining activity. As a result of this subsidence and subsequent soil erosion, trees cannot live in these areas.

MODIS Data
The MODIS vegetation index product, MOD13Q1, with a 16-day composite and a spatial resolution of 250 m, was selected to provide temporally frequent coarse resolution data. It contains both the NDVI product and reflectance layers of near infrared (NIR) and red that will be used in the following analysis. All MOD13Q1 scenes of h26v04 and h26v05 from 2000 to 2011 were obtained from the National Aeronautics and Space Administration (NASA) Earth Observing System Data and Information System (http://earthdata.nasa.gov). To meet the STARFM input requirements, MOD13Q1 data was subset to the extent of the study area, re-projected to the same Universal Transverse Mercator (UTM) projection as Landsat and resampled to a resolution of 30 m using the MODIS re-projection tool (MRT). A nearest neighbor resample approach was used to ensure the least distortion of the original NDVI value.
The MOD13Q1 vegetation product contains a summary quality layer called "pixel reliability", holding information on the overall pixel quality. In this study, only pixels with flag values equal to 0 (good data) were used, whereas pixels with other flag values were masked out. A detailed explanation of MOD13Q1 pixel reliability can be found at http://www.ctahr.hawaii.edu/grem/mod13ug/ sect0005.html. As MODIS data were not collected in the first two months of the year 2000, average values of the same period between 2001 and 2011 were used to form a completed time series that was needed for the trend analysis [41].
A simple five-point temporal filter was applied to the MOD13Q1 NDVI time series to fill gaps and missing values from bad quality pixels in the original time series [41][42][43]. For a missing NDVI value at position X in the time series, the average of NDVI values at X − 1 and X + 1 is given. If only one of the two values is available, then this value will be used. If both the values are missing, NDVI values at X − 2 and X + 2 will be included in the filtering process, following the same rules as for X − 1 and X + 1.

Landsat Data
Landsat TM/ETM+ images L1T (Level 1 Terrain Corrected) products of WRS (Worldwide Reference System)-2, Path 125/Row 32 were used. A total number of 32 Landsat scenes were available (cloud cover <10% and SLC (Scan Line Corrector) -on) from 2000 to 2011 and all of them were downloaded using the United States Geological Survey (USGS) GLOVIS (Global Visualization Viewer) portal (http://glovis.usgs.gov). These Landsat images were found to be spatially consistent at a high level in the study area, thus no additional registration work was done. Atmosphere correction of Landsat images was performed using the LEDAPS (Landsat Ecosystem Disturbance Adaptive Processing System) tool [44], in which column water vapor is taken from NOAA (National Oceanic and Atmospheric Administration) National Centers for Environmental Prediction (NCEP) reanalysis data, ozone concentrations are derived from Total Ozone Mapping Spectrometer (TOMS) data, and aerosol optical thickness (AOT) is directly extracted from the imagery using the Dark Dense Vegetation (DDV) method [45]. As the LEDAPS employs the same radiative transfer model for MODIS atmospheric correction, the difference between Landsat and MODIS reflectance data was minimized. The NIR and red reflectance bands in Landsat were then extracted and subset to the extent of the study area. Under ideal conditions, atmospheric-corrected Landsat observations should be identical to the MODIS surface reflectance products, which is one of the assumptions of the STARFM algorithm. However, systematic deviation between the two datasets will always exist [31,33] that could affect the performance of STARFM. To minimize this impact, Landsat NIR and red reflectance bands were normalized to the corresponding MODIS images by subtracting the mean difference value between them.

STARFM
Neglecting geolocation errors and differences in atmospheric correction, a coarse-resolution MODIS pixel can be aggregated from the corresponding fine-resolution Landsat pixels. For a homogenous pixel at the MODIS resolution, Landsat surface reflectance can be obtained by adding an error value (caused by differing bandwidth and solar geometry) to the MODIS surface reflectance [21]. Suppose the land cover type and system errors do not change over time; we will be able to predict the Landsat surface reflectance on a specific day (t n ) by blending one or more pairs of Landsat and MODIS surface reflectance images on a base day (t 0 ) and the MODIS surface reflectance image acquired on day t n [21]. However, the relationships between MODIS and Landsat reflectance are complicated by the fact that MODIS observations may be mixed pixels, and land cover type and Bidirectional Reflectance Distribution Function (BRDF) may change over time [21]. To solve this problem, STARFM employed a moving window to incorporate additional information from neighboring Landsat and MODIS image pixels using a weighting function defined by a combination of spectral, temporal and spatial differences among the given Landsat and MODIS scenes [21]. Surface reflectance time series can be obtained at the same spatial resolution as Landsat and the same temporal resolution as MODIS using the STARFM algorithm. The synthetic work was performed using the STARFM tool acquired from the Landsat Ecosystem Disturbance Adaptive Processing System (http://ledapsweb.nascom.nasa.gov/).

Two Schemes of using STRFM
Two schemes were proposed to generate a synthetic time series of 30-m resolution NDVI images using STARFM ( Figure 3). Scheme 1: The NIR and red reflectance bands of Landsat and MODIS images on day t 0 and MODIS bands (NIR and red) on day t n were used as inputs for STARFM to generate synthetic bands (NIR and red) on day t n . Then, the NDVI image on day t n was calculated using the synthetic NIR and red bands generated above. Scheme 2: the Landsat NDVI image on day t 0 was first derived from NIR and red bands, which was then used directly as an input for STARFM, together with MODIS (MOD13Q1) NDVI images on day t 0 and t n , to generate the synthetic NDVI image on day t n . Both of these two schemes were used to generate synthetic NDVI images, and the accuracy of the results were evaluated by comparing with Landsat NDVI on day t n . A random sample (25%) of the vegetation pixels in both independent Landsat and synthetic NDVI images were compared on a per-pixel basis by means of correlation. The scheme with higher accuracy was chosen to generate synthetic NDVI time series from 2000 to 2011 and was further used for the following trend analysis. Table 1 shows the base pair images used in STARFM to generate the synthetic NDVI time series from 2000 to 2011. Although 32 Landsat scenes were available in this study area, nearly half of them were obtained in the non-growing season, and the temporal distribution of these scenes was very unbalanced. Consequently, one base pair of Landsat and MODIS images were selected for each year. For those years where no Landsat images were available, base pair images from the adjacent year were used.

Accuracy Assessment of STARFM Synthetic Data
Available Landsat images in 2000 were used for accuracy assessment of the two schemes above, as there were five cloud-free Landsat scenes covering the whole growing season (6 May  To evaluate the accuracy of synthetic NDVI time series, the mean and standard deviation values of synthetic dense NDVI time series images and NDVI max time series images were calculated and compared with the corresponding MODIS time series.

Trend Analysis and Evaluation
Trend analysis of NDVI time series was carried out in two stages. First, all of the synthetic NDVI images were included in the analysis to ensure that the temporal properties of the synthetic time series data were fully utilized. Theil-Sen (TS) slope trend analysis was used to quantify the temporal trends in the NDVI time series. This is a non-parametric statistical method calculating the median of the slopes between all pair-wise combinations over time [46,47]. The TS slope can tolerate up to approximately 30% noisy input observations without influencing the trend estimate [48] and is robust against seasonality, non-normality and serial dependence [49,50]. The significance of the NDVI time series trend was examined by the non-parametric seasonal Mann-Kendall (SMK) test, which is a completely rank-based test that is robust against non-normality, missing values, seasonality and, if corrected, can be robust against serial dependence [49,51,52]. The SMK test produces Z scores to measure the significance of a monotonic trend over time, with a Z value >1.96 representing a significant increasing trend and a Z value <−1.96 indicating a significant decreasing trend (p < 0.05).
In this study, SMK Z scores were obtained using the method proposed by Hirsch et al. [49], considering serial dependence in the seasonal time series. This method was recommended by de Beurs et al. [51] for remote sensing time series analysis and has been used for vegetation trend analysis with dense time series satellite data [53][54][55].
In the second stage, the annual maximum NDVI (NDVI max ) values for every pixel in each year between 2000 and 2011 were extracted from both synthetic and MODIS time series data and used for trend analysis. NDVI max values represent the highest vegetation productivity in one year, thus avoiding seasonal variation and minimizing the errors introduced by STARFM. In this stage, temporal trends in the NDVI max time series were also examined by the TS slope. The Mann-Kendall (MK) test [56,57], which does not consider seasonal variation, was used to calculate the significance of the trends.
Both the NDVI and NDVI max trends derived from STARFM synthetic time series were evaluated by comparing to the corresponding trends derived from the MODIS time series.   Figure 4A shows the results of scheme 1 using band (NIR and red) reflectance images on 13 October 2000 as inputs for STARFM, while Figure 4B shows the results of scheme 2 using NDVI images on 13 October 2000 directly as an input for STARFM. For the four validation images in this study area, the results of scheme 2 showed a higher accuracy, with scatterplots closer to the 1:1 line and R 2 > 0.70, in contrast to those of scheme 1, in which 0.56 < R 2 < 0.70. The figures for different vegetation types in Table 2 give similar results, and scheme 2 consistently generates higher R 2 values and lower mean absolute difference values than scheme 1. Consequently, scheme 2 was selected to generate synthetic NDVI time series from 2000 to 2011.

Synthetic NDVI Time Series
To  Table 2 for R 2 values and mean absolute difference for different vegetation types. Therefore, with base pair images from the adjacent year, STARFM could still maintain the accuracy of synthetic NDVI images. Synthetic and MODIS NDVI images were compared in terms of the mean and standard deviation of dense time series images and NDVI max time series images during the period of 2000-2011 ( Figure 5). The mean values in the two time series showed a high level of consistency, while synthetic NDVI images show slightly higher standard deviation values than MODIS images. It is likely that synthetic images maintained the high spatial resolution characteristics from Landsat images and captured, at the same time, NDVI changes in the MODIS time series. The NDVI max images using the same base Landsat image (Table 1) showed similar standard deviation values ( Figure 5B), possibly since the spatial patterns of the synthetic images were determined by the base Landsat image [21].

Comparison of Trends from Synthetic and MODIS Time Series
The trends from NDVI dense time series and annual NDVI max time series during the period of 2000-2011 were obtained from both synthetic and MODIS (MOD13Q1) time series at a spatial resolution of 30 m and 250 m, respectively. Using MODIS as the reference time series, comparisons between the two results were carried out, to evaluate the accuracy of the NDVI trends derived from the synthetic time series. Figure 6 shows the spatial distribution of NDVI trends derived from both synthetic and MODIS dense time series for the entire study area during the period of 2000-2011. Figure 6A,B shows the magnitude of NDVI trends over time (NDVI slopes shown in NDVI units per year), while Figure 6C,D shows the statistical significance of the trends, according to the SMK test (Z scores). The results of the synthetic time series shown in Figure 6A,C exhibit a similar overall pattern with the results of the MODIS time series presented in Figure 6B,D. With a spatial resolution of 30 m, however, the results of the synthetic time series provide more detailed information than the results of the MODIS time series, with a spatial resolution of 250 m. More areas with a negative trend were found in the results of the synthetic time series than in the results of the MODIS time series, especially in the mined area, although positive trends dominated for both time series.  Figure 7 shows the scatterplots of the SMK Z scores versus NDVI slopes for the vegetated pixels in the study area, indicating the relationship between the magnitude of NDVI trends and the statistical significance of the trends. The total number of pixels in the NDVI trends derived from the synthetic time series is nearly 70-times larger than in those derived from the MODIS time series, which equals the ratio of per pixel area between MODIS and Landsat images. In both scatterplots, a small proportion of pixels have a positive slope, but a negative Z value (between −1.96 and zero in the second quadrant). The discrepancy between TS and SMK was also reported by Fensholt et al. [54], who used a prewhitening method [58] to correct the serial correlation within the seasonal time series. However, this inconsistency in the results of the two statistical methods could be tolerated, because none of the positive TS slopes in the second quadrant passed the SMK significance test. This explains the reason for the more positive trends in the NDVI slope maps ( Figure 6A,B) than in the corresponding significance maps (Figure 6C,D). The majority of the high frequency pixels in Figure 7A show a similar pattern with those in Figure 7B in terms of the relationship between Z scores and TS slopes, while the NDVI slopes in Figure 7A extend to wider ranges. This may be caused by the higher spatial resolution of the synthetic images, which allows them to capture NDVI changes at the scale smaller than the spatial resolution of the MODIS images. The frequency distributions of NDVI slopes in Figure 6A,B were further studied and are shown in Figure 8. The mean values of the NDVI slopes of both NDVI trends were nearly the same, with 0.0036 for those derived from synthetic time series and 0.0031 for those derived from MODIS time series. However, the standard deviation value was much higher for NDVI slopes derived from synthetic time series (0.0057) than those derived from MODIS time series (0.0020). Due to the huge difference of the standard deviation values and the total pixel numbers between these two distributions, both the means and standard deviation values of them were found to be significantly (p < 0.001) different from each other according to the t-test (two-side) and F-test, respectively. In both histograms, pixels characterized with slope values between 0.002-0.004 account for the largest proportion. However, the percentages of these pixels differed substantially, with 27.6% of the total vegetated pixels in the NDVI slopes concentrated between 0.002-0.004 derived from the synthetic time series, while 63.1% in that derived from the MODIS time series. Possible reasons for these differences are discussed in Section 4.2.  Table 3 reflects similar spatial patterns in the two NDVI trends, showing a cross-tabulation with the vegetated surface area characterized by positive and negative trends, according to the SMK test. The p-value obtained from a chi-square test of the cross-tabulation table was less than 0.001, meaning that there was a statistically significant relationship between the two trends [59]. The contingency coefficient, which is recognized as an objective measure for the degree of dependence of the classifications in the cross-tabulation table, was 0.56. For a 4 × 4 matrix, the maximum contingency coefficient is 0.86. Consequently, the NDVI trends derived from synthetic dense time series showed a high level of consistency with the NDVI trends derived from MODIS time series, in terms of the significance of the NDVI trends obtained by the SMK test. Table 3. Cross-tabulation considering the significance of NDVI trends derived from the synthetic and MODIS dense time series (2000-2011). Chi-square test of the cross-tabulation showed a significant (p < 0.001) relationship between the two trend maps with a contingency coefficient of 0.56 out of 0.86 (values in km 2 ).  Figure 9 shows the spatial distribution of the NDVI max trends derived from both synthetic and MODIS NDVI max time series during the period of 2000-2011. Figure 9A,B shows the magnitude of NDVI max trends over time (NDVI max slopes in NDVI units per year), while Figure 9C,D shows the statistical significance of the trends, according to the MK test (Z scores). The patterns between the two NDVI max trends were visually consistent at a high level.  Figure 10 shows the scatterplots of the MK Z scores versus NDVI max slopes for the vegetated pixels in the study area, and similarities in the patterns of the relationship between Z scores and NDVI max slopes were evident. All of the pixels are located in the first and third quadrants, showing the consistent results of the TS slope estimator and the MK significance test. The frequency distributions of NDVI max slopes in Figure 9A,B are shown in Figure 11. The mean values of NDVI max slopes of both trends were nearly the same, with 0.0009 for those derived from synthetic time series and 0.0007 for those derived from MODIS time series. The standard deviation was higher for NDVI max slopes derived from synthetic time series (0.0090) than for those derived from MODIS time series (0.0065). Like the results of Figure 8, both the means and standard deviation of these two distributions were significantly (p < 0.001) different according to the t-test (two-side) and F-test, respectively. The possible reasons for these differences are discussed in Section 4.2.  The cross-tabulation of the results of the NDVI max trends is shown in Table 4, illustrating the vegetated surface area characterized by positive and negative trends, according to the MK test. A statistically significant relationship was found between the two NDVI max trends, with a p-value less than 0.001 and a contingency coefficient of 0.57 out of 0.86 obtained from the chi-square test of the cross-tabulation table. Therefore, the NDVI max trends derived from synthetic time series were highly consistent with the NDVI max trends derived from the MODIS time series, in terms of the significance of NDVI max trends obtained by the MK test. Table 4. Cross-tabulation considering the significance of NDVI max trends derived from the synthetic and MODIS time series (2000-2011). The chi-square test of the cross-tabulation showed a significant (p < 0.001) relationship between the two trend maps with a contingency coefficient of 0.57 out of 0.86 (values in km 2 ).

Accuracy of the Synthetic NDVI Images
The accuracy of STARFM synthetic NDVI images derived from scheme 2, in which NDVI images were used directly as input, was higher than that derived from scheme 1, in which band reflectance images were used as input ( Figure 4A,B). This may be because the nonlinear relationships between NDVI and reflectance bands amplified the errors in the synthetic NIR and red reflectance images when transforming into NDVI images, which was also reported in the results of Walker et al. [32]. However, the decrease of synthetic accuracy in NDVI was seldom seen in the results of Bhandari et al. [33] in terms of the R 2 values, indicating that the performance of the two schemes may differ between separate study areas. The huge numbers of ravines and gullies in the study area increase the variation of BRDF (Bidirectional Reflectance Distribution Function) and the heterogeneity of the land surface, which may explain the lower accuracy of synthetic NDVI images ( Figure 4, Table 2) compared to other studies, where the landscapes were more flat and homogeneous [30,32,33]. Bhandari et al. [60] examined the impact of solar and viewing geometry on the MOD13Q1 product in an Australian forest and found that there existed a higher amplitude and a phase shift in the MOD13Q1 NDVI time series compared to NDVI derived from the BRDF adjusted MODIS product (MCD43A4). Hence, using the MOD13Q1 product as input for STARFM may introduce errors in the synthetic NDVI images, but the MOD13Q1 product still has an important advantage over the MCD43A4, as the spatial resolution of MOD13Q1 is the highest among MODIS products. In addition, performing topographic correction on the Landsat and MODIS images may improve the results of STARFM outputs, especially in mountainous regions.
The spectral difference between pixels in base pair Landsat and MODIS images, regarded as a measure of the homogeneity of the MODIS pixel, is one of the three factors that decide the weight function in STARFM [21]. Thus, spectral differences caused by other reasons, such as different sensors and spectral response function, will lower the accuracy of synthetic images or increase the differences of the mean values between synthetic and MODIS images [33]. Therefore, the normalization of Landsat reflectance images could minimize these kinds of errors to obtain a consistent synthetic time series with MODIS.
In the study of Bhandari et al. [33] and Schmidt et al. [35], the number of base Landsat images used to generate the long-term synthetic time series were much larger than the number used in this study, with 38 scenes in five years and 90 scenes in 7.5 years, respectively. We only found 32 Landsat scenes in 12 years in this study area and, due to the unbalanced dates of observation and many scenes from the non-growing period, only six of them were finally used to generate the synthetic time series in our study. When using different base Landsat images with STARFM, Walker et al. [32] reported minor differences between the synthetic images on the same date. Therefore, despite the small amount of base Landsat images in our study, the accuracy of our synthetic time series were not significantly decreased. In addition, all of the six base Landsat images used in this study were obtained within the same season (between 21 August and 14 October), ensuring the same accuracy level of synthetic images at the same date in each year, which is important for the SMK test.

Potentials of the Trends Derived from Synthetic Time Series
At a spatial resolution of 250 m, NDVI trends derived from MODIS time series only captured the main structures of vegetation changes, while the fine-scale patterns of negative and positive trends were neutralized or even eliminated ( Figures 6B,D and 9B,D). In contrast, at a spatial resolution of 30 m, synthetic time series could capture finer patterns of positive and negative trends at the sub-MODIS pixel scale in Figure 6A,C and Figure 9A,C. These results were supported by the study of Stellmes et al. [8], who examined the effects of spatial resolution on NDVI trends by spatially degraded Landsat images. Positive and negative slopes were neutralized in the trend map derived from MODIS time series, due to the large pixel size, while these two different slopes could be separated in the synthetic time series, thus reducing the proportion of pixels with a moderate slope value and enlarging the ranges (Figures 8 and 11). Therefore, synthetic time series captured vegetation change information that could not be captured by MODIS time series, making the distributions of NDVI slopes statistically different (p < 0.001) between the two time series. Consequently, STARFM-generated synthetic NDVI time series have the potential to enhance the ability to capture finer scale patterns of vegetation changes.
According to the statistical significance test (SMK and MK), the trends calculated from synthetic dense NDVI time series and NDVI max time series were different. Specifically, 68% of the vegetated areas were characterized by positive trend (41% being significant and 27% non-significant) in the synthetic dense time series, while the figure for MODIS dense time series was even higher, at 87% (56% being significant and 31% non-significant). By contrast, in the NDVI max time series, areas characterized by positive and negative trends were roughly at the same level, with 53% experiencing a positive trend (9% being significant and 43% non-significant) in synthetic series and 52% for MODIS series (5% being significant and 48% non-significant). The magnitudes of both positive and negative trends in NDVI max time series are higher than those in dense time series. These indicate that different vegetation trends were signified with different time series data around the Datong Coalfield. For example, an area may experience a positive trend in NDVI time series for the spring, autumn and winter seasons, but a negative one for the summer season, which may result in an increasing overall trend by examining dense NDVI time series with SMK test and TS slope for this area, but a decreasing trend by NDVI max time series with MK test and TS slope. Lasanta and Vicente-Serrano [9] reported different trend directions between August and March time series in northeast Spain, which supports the above assumption. This may explain the different trend directions between the dense NDVI time series and NDVI max time series described above. Consequently, it is worthy to perform trend analysis for NDVI time series of particular seasons or months, which may not be reflected in dense time series or NDVI max time series. STARFM also provides a feasible way to do such a job at a spatial resolution of 30 m, which is especially valuable in regions where the number of available Landsat scenes is very limited. However, field works need to be conducted to validate the trends obtained from STAFRM-generated time series data and to decide which time series is the most appropriate one on quantifying mining effects on vegetation.
STARFM usually generates poorer synthetic images in heterogeneous areas than in homogeneous areas [21], thus a lower consistency level or even totally different directions of the trends may be obtained between synthetic and MODIS time series in heterogeneous areas. The potential use of generating high spatial and high temporal resolution synthetic NDVI time series is to better quantify the effects of mining activities on vegetation by comparing vegetation trends between mined and non-mined areas. Therefore, the time series data generated in those areas should be accurate enough (i.e., consistent with MODIS time series) to obtain reliable results. Areas characterized by significantly different trends between synthetic and MODIS time series should not be selected to conduct such comparison work. An example would be the central part of the mined area, where the decreasing NDVI trend was obvious in the synthetic dense time series in Figure 6A,C, while an opposite NDVI trend was seen in the MODIS dense time series in Figure 6B,D. However, the decreasing NDVI max trend of this area was consistent in both the synthetic and MODIS time series (Figure 9), indicating that the synthetic NDVI max time series data were able to obtain consistent results in this area. This may be because errors in the synthetic images were reduced by calculating the annual maximum NDVI images.
Synthetic images generated by the same base Landsat scene were highly inter-correlated (Table 1, Figure 5), which was also reported by Schmidt et al. [35]. Therefore, the scarcity of Landsat scenes in the study area added the strong serial dependence of the synthetic time series, as they were produced by STARFM using base pair images from the adjacent year. Thus, the statistical methods that are robust against serial dependence used in this study were critical to obtain reliable NDVI trends from the synthetic dense time series.

Conclusions
Extensive coal mining activities pose a threat to the local environmental conditions of the Loess Plateau, where the ecosystems are very fragile. To accurately quantify the vegetation changes caused by coal mining activates, this study proposed a method to generate high spatial and high temporal time series NDVI data using the Spatial and Temporal Adaptive Reflectance Fusion Model (STARFM). The main conclusions can be summarized as follows: (1) Using Landsat and MODIS NDVI images directly as input for the STARFM algorithm generated more accurate synthetic NDVI images than using band reflectance images as input for the STARFM algorithm around the Datong Coalfield. (2) It is feasible to use one base pair image to generate synthetic time series for a year and use base pair images from the adjacent years when no Landsat scenes are available in a certain year. (3) NDVI trends from synthetic dense time series (2000-2011) were found to be highly consistent with that from MODIS time series around the Datong Coalfield, while the consistency level was even higher for annual maximum NDVI (NDVI max ) trends from synthetic and MODIS time series. (4) Despite the fact that the overall patterns of NDVI trends were maintained, trends from MODIS time series at a spatial resolution of 250 m neutralize fine-scale positive and negative vegetation changes. Synthetic time series generated via STARFM are likely to be more capable in capturing such fine-scale vegetation changes. (5) It is suggested that when using the synthetic time series to quantify the effects of mining activities on vegetation, the test areas should be selected with caution, as the NDVI trends obtained from synthetic high spatial resolution time series and MODIS time series may be significantly different in some heterogeneous areas.