Detecting Historical Vegetation Changes in the Dunhuang Oasis Protected Area Using Landsat Images

: Given its proximity to an artiﬁcial oasis, the Donghu Nature Reserve in the Dunhuang Oasis has faced environmental pressure and vegetation disturbances in recent decades. Satellite vegetation indices (VIs) can be used to detect such changes in vegetation if the satellite images are calibrated to surface reﬂectance (SR) values. The aim of this study was to select a suitable VI based on the Landsat Climate Data Record (CDR) products and the absolute radiation-corrected results of Landsat L1T images to detect the spatio-temporal changes in vegetation for the Donghu Reserve during 1986–2015. The results showed that the VI difference ( ∆ VI) images effectively reduced the changes in the source images. Compared with the other VIs, the soil-adjusted vegetation index (SAVI) displayed greater robustness to atmospheric effects in the two types of SR images and was more responsive to vegetation changes caused by human factors. From 1986 to 2015, the positive changes in vegetation dominated the overall change trend, with changes in vegetation in the reserve decreasing during 1990–1995, increasing until 2005–2010, and then decreasing again. The vegetation changes were mainly distributed at the edge of the artiﬁcial oasis outside the reserve. The detected changes in vegetation in the reserve highlight the increased human pressure on the reserve.


Introduction
Dunhuang was an important stop on the ancient Silk Road.Today, it is a famous cultural heritage city.This city is located in Northwestern China in the Western Hexi Corridor in Gansu Province; it lies within the triangle formed by Gansu, Qinghai, and Xinjiang Provinces [1][2][3][4][5][6][7][8][9].In recent decades, the water-based geological and environmental problems of the Dunhuang basin have worsened, as indicated by the decline of its wetlands and the degradation of its vegetation.The Dunhuang Xihu, Nanhu, Beihu, and Donghu Nature Reserves were established to maintain a natural "green barrier" to protect the ecological diversity and environment of the desert oasis [3,4,10].The Donghu Nature Reserve is located to the east of Dunhuang and near the artificial oasis; it connects between the Guazhou Oasis and the Dunhuang Oasis (Figure 1) and is of considerable importance to the stability of the two oases.The increasing environmental pressure on the Donghu Nature Reserve is primarily due to uncontrolled human exploitation of the area and its water resources.In the 1960s and 1970s, the upper reaches of the Shule River and the Dang River were dammed.This intervention caused sections of the rivers to be cut off, reduced the area of wetlands, and caused the decline and die-off of natural vegetation in the reserve [5][6][7][8][9].In addition, continued increases in the population, the area of the artificial oases [11], and the large-scale planting of cotton, grapes, and other crops that consume large amounts of water in Dunhuang caused a surge in agricultural irrigation.Much of the groundwater has been exploited, leading to a continuing decline in groundwater levels.The problem of water shortages in the nature reserve has become increasingly prominent [3,4].Protected areas are important for the sustainable economic development of the region [5,8].The environmental problems of the protected areas in Dunhuang have attracted widespread attention.However, most of the research on the protected areas has focused on the Xihu Reserve [5][6][7][8] and Nanhu Reserve [9], with few studies addressing the Donghu Reserve.Detection of historical changes in the vegetation of the Donghu Nature Reserve over the last 30 years could help the government understand the historical landscape of the reserve.Understanding these changes in vegetation would facilitate identifying the drivers of recent vegetation changes in the Donghu Reserve and provide a reference for research into the vegetation of other reserves.
Vegetation is the primary component of terrestrial ecosystems, and it is a natural medium that connects soil, air, and moisture.It also represents the general status of the environment in a region and acts as an indicator that is used in global change research [12,13].Vegetation indices (VIs) are important remote sensing (RS) parameters because they reflect the status of vegetation through the digital combination of information from diverse spectrum bands, especially the visible and nearinfrared bands.VIs reflect the collective status of chlorophyll, leaf area, coverage, and canopy structure [14,15].Repeated VI images of the same area can be used to measure variations in the biophysical characteristics of vegetation, including vegetation coverage.The difference between a later image and an earlier image is called the VI difference (ΔVI) [16].In a ΔVI image, the unchanged pixels surround the mean value of 0, whereas the changed pixels cluster at the positive and negative ends of the distribution.
Many external factors must be considered when analyzing images of the same area taken at different times; these factors include the soil background, atmospheric conditions, topography, illumination, observation angles, and sensor calibration, all of which can affect the Vis' values [17,18].The increasing environmental pressure on the Donghu Nature Reserve is primarily due to uncontrolled human exploitation of the area and its water resources.In the 1960s and 1970s, the upper reaches of the Shule River and the Dang River were dammed.This intervention caused sections of the rivers to be cut off, reduced the area of wetlands, and caused the decline and die-off of natural vegetation in the reserve [5][6][7][8][9].In addition, continued increases in the population, the area of the artificial oases [11], and the large-scale planting of cotton, grapes, and other crops that consume large amounts of water in Dunhuang caused a surge in agricultural irrigation.Much of the groundwater has been exploited, leading to a continuing decline in groundwater levels.The problem of water shortages in the nature reserve has become increasingly prominent [3,4].Protected areas are important for the sustainable economic development of the region [5,8].The environmental problems of the protected areas in Dunhuang have attracted widespread attention.However, most of the research on the protected areas has focused on the Xihu Reserve [5][6][7][8] and Nanhu Reserve [9], with few studies addressing the Donghu Reserve.Detection of historical changes in the vegetation of the Donghu Nature Reserve over the last 30 years could help the government understand the historical landscape of the reserve.Understanding these changes in vegetation would facilitate identifying the drivers of recent vegetation changes in the Donghu Reserve and provide a reference for research into the vegetation of other reserves.
Vegetation is the primary component of terrestrial ecosystems, and it is a natural medium that connects soil, air, and moisture.It also represents the general status of the environment in a region and acts as an indicator that is used in global change research [12,13].Vegetation indices (VIs) are important remote sensing (RS) parameters because they reflect the status of vegetation through the digital combination of information from diverse spectrum bands, especially the visible and near-infrared bands.VIs reflect the collective status of chlorophyll, leaf area, coverage, and canopy structure [14,15].Repeated VI images of the same area can be used to measure variations in the bio-physical characteristics of vegetation, including vegetation coverage.The difference between a later image and an earlier image is called the VI difference (∆VI) [16].In a ∆VI image, the unchanged pixels surround the mean value of 0, whereas the changed pixels cluster at the positive and negative ends of the distribution.
Many external factors must be considered when analyzing images of the same area taken at different times; these factors include the soil background, atmospheric conditions, topography, illumination, observation angles, and sensor calibration, all of which can affect the Vis' values [17,18].These factors often contribute substantial amounts of noise, which can affect the results of applying the VI values.To accurately evaluate changes in a VI, the images must be calibrated to surface reflectance (SR) data to eliminate or reduce the influence of atmospheric conditions [19,20].The methods used to invert the true reflectivity of features include both relative and absolute radiation corrections.The FLAASH model is based on the MODTRAN5 radiation transmission model and is commonly used to perform absolute radiation corrections.This model obtains SR values with relatively high accuracy [21].Landsat's advanced CDR products provides SR images from the TM sensor [22].The provisional Landsat 8 Surface Reflectance Code (LaSRC), which employs Moderate Resolution Imaging Spectroradiometer (MODIS) data, is distinctly different from the algorithm used by the USGS to process TM L1T products to obtain SR values [23,24].
The overall goal of the present study was to detect spatio-temporal changes in vegetation in the Donghu Reserve using seven Landsat datasets from different time points (1986,1990,1995,2000,2005,2010,2015).The specific objectives were: (1) explore the impact of the differences between the VIs produced using the CDR products and the results of absolute radiation corrections of the Landsat L1T images on the inferred vegetation changes in the reserve, (2) determine a suitable VI for assessing the vegetation changes within the study area, and (3) analyze the historical changes in the vegetation of the reserve from a spatio-temporal perspective.The present study provides a basis for the study of vegetation changes in other reserves.

Study Area
The Donghu Nature Reserve (40 • 7 51"-40 • 28 58" N, 94 • 47 45"-95 • 19 41" E) is located along the Dang River Plain in Dunhuang on the eastern side of an artificial oasis.It covers Yitang Lake, Xindian Lake, and Daquanwan and has an area of approximately 2.7 × 10 4 km 2 .The average elevation of the reserve is 1110 m, and the average annual precipitation and evapotranspiration are 39.9 mm and 2490 mm, respectively.The reserve has a warm temperate arid climate.The abundances of water and plants represent barometers of the ecological environment of the desert oasis in Dunhuang, which has been well protected by the local government and relevant departments for many years.However, because of its proximity to the artificial oasis, the reserve has been affected by drought and agricultural production activities.Serious desertification has been observed in the area surrounding the reserve, reeds and other wetland vegetation have declined, and the plant community has experienced serious degradation [25].In addition, the Yitanghu wetland is rich in nitrate and thick shrubs, and many springs are present in the surrounding area; thus, this area is subject to frequent human activities, such as mining, grazing, digging, and hunting.These activities have exacerbated the degradation of the vegetation in the reserve [26].The boundaries of the Donghu Reserve were extracted based on the forest resources distribution map of Dunhuang (1:250,000) for 2004.In this study, because the vegetation was the main object of analysis, the mountains in the reserve were masked out, and the area within a 2-km buffer surrounding the reserve (shown by the red line in Figure 1) was the main study area.

Data
The Landsat images used in this study, which cover the period from 1986 to 2015 and are located at path 137 and row 32, were downloaded from the USGS.The scenes, which were collected by the Landsat 5 TM and Landsat 8 OLI satellites, have a resolution of 30 m and are terrain-corrected Level 1T products and the corresponding Landsat CDR-SR products.As shown in Table 1, the time period covered for each year was July to August [11].This period corresponds to the best phenological period for the vegetation.The use of a consistent period avoided differences caused by different seasons and allowed for the differentiation of different features in the images.

Method
The workflow, shown in Figure 2, primarily includes data preprocessing, evaluation of the FLAASH atmospheric correction, the calculation of the VIs and ∆VIs, the extraction of the vegetation changes based on threshold criteria and the production of maps of the spatio-temporal vegetation changes.These steps are described in detail below.

Method
The workflow, shown in Figure 2, primarily includes data preprocessing, evaluation of the FLAASH atmospheric correction, the calculation of the VIs and ΔVIs, the extraction of the vegetation changes based on threshold criteria and the production of maps of the spatio-temporal vegetation changes.These steps are described in detail below.

Data Preprocessing
To ensure that pairs of images occupying the same position had consistent coordinates, image registration of the different phases was performed when possible to eliminate "pseudo-changes" caused by geometric errors [27].A Landsat CDR-SR image obtained on 8 August 2015, was selected as the reference image; as required, this image largely reflected clear atmospheric conditions.Although this image did contain cloud-covered areas, the pixels within the study area were not affected [28].Based on the "Map_Registration_Select GCPs: Image to Image" module in software ENVI v5.1 (Exelis Visual Information Solutions, Boulder, CO, USA), 42 uniformly-distributed points were selected for image registration.To ensure the accuracy of the results, surface pixels that were not prone to warpage and spectral deformation, such as road intersections and inflection points, river confluences, and buildings that did not change over time, were selected as control points.Using the second-order polynomial correction, the registration errors were found to be equal to or less than 0.5 pixels.
The digital number (DN) of each image was transformed into radiance by radiometric calibration.The following conversion formulas were applied to the TM and OLI data, respectively: where QCAL represents the DN of a pixel;

Data Preprocessing
To ensure that pairs of images occupying the same position had consistent coordinates, image registration of the different phases was performed when possible to eliminate "pseudo-changes" caused by geometric errors [27].A Landsat CDR-SR image obtained on 8 August 2015, was selected as the reference image; as required, this image largely reflected clear atmospheric conditions.Although this image did contain cloud-covered areas, the pixels within the study area were not affected [28].Based on the "Map_Registration_Select GCPs: Image to Image" module in software ENVI v5.1 (Exelis Visual Information Solutions, Boulder, CO, USA), 42 uniformly-distributed points were selected for image registration.To ensure the accuracy of the results, surface pixels that were not prone to warpage and spectral deformation, such as road intersections and inflection points, river confluences, and buildings that did not change over time, were selected as control points.Using the second-order polynomial correction, the registration errors were found to be equal to or less than 0.5 pixels.
The digital number (DN) of each image was transformed into radiance by radiometric calibration.The following conversion formulas were applied to the TM and OLI data, respectively: where QCAL represents the DN of a pixel; QCAL max and QCAL min represent the maximum and minimum DN values of the image, respectively; and L max and L min represent the maximum and minimum radiance measured by the sensor, respectively.These values were determined from the metadata of each image.

FLAASH Atmospheric Correction
FLAASH is an atmospheric correction model that is based on the MODTRAN5 radiation transmission model and produces high-precision results.The atmospheric properties were estimated via feature points on the spectrum of the image pixels using the model instead of using atmospheric data collected at the same time that the images were taken.This process effectively removed scattering by vapor and aerosols.The pixel-based corrections were also able to rectify the effects of cross radiance caused by the proximity between the target pixels and adjacent pixels [29].In addition, the model performed effective spectral smoothing of the noise caused by artificial suppression, and accurate parameters of the true physical model, such as reflectivity, emissivity, and surface temperature, were obtained.The atmospheric corrections used in the FLAASH model are based on the standard planar Lambertian (or the approximate planar Lambertian) in the range of the solar spectrum in addition to the thermal radiation [21,29].In this study, seven images with absolute atmospheric corrections were produced using the FLAASH module embedded in ENVI 5.1.

Vegetation Indices
The following four main VIs were tested to determine the applicability of the CDR-SR and the SR calculated with the FLAASH model (Fa-SR) images: the Normalized Difference Vegetation Index (NDVI), the Generalized Difference Vegetation Index (GDVI), the Global Environment Monitoring Index (GEMI), and the Soil-Adjusted Vegetation Index (SAVI).These VIs employ the spectral information contained in the red (ρ R ) and near-infrared (ρ NIR ) bands, which facilitates the delineation of vegetation compared with the delineation obtained using other bands.The NDVI and GDVI are ratio indices and derivatives of the Ratio Vegetation Index (RVI) [15].The SAVI and GEMI use fixed constants instead of the band ratio model.The NDVI is the most widely used VI [30].The SAVI can reduce the impact of the soil background in areas of sparse vegetation [31], whereas the GDVI has higher sensitivity in areas of sparse vegetation [14].The GEMI performs well in detecting vegetation in arid areas [32,33].
Since the differences between the CDR-SR images and the Fa-SR images could potentially affect the VI images, a correlation analysis and a linear regression analysis of the VI CDR images (which were produced using the CDR-SR images) and the VI Fa images (which were produced using the Fa-SR images) were performed to assess the sizes of these differences.

Detecting Changes in the Vegetation
The VI CDR images were ordered in chronological order, and the ∆VI CDR images were produced by taking the subsequent image (t 2 ) minus the previous image (t 1 ) in pairs of consecutive images [22].Here, each ∆VI CDR included the four ∆VIs produced from the four Vis' Equation (3): ∆GDVI CDR , ∆GEMI CDR , ∆NDVI CDR , and ∆SAVI CDR .Similarly, each ∆VI Fa image included the four ∆VIs produced from the four VIs' Equation (4): ∆GDVI Fa , ∆GEMI Fa , ∆NDVI Fa , and ∆SAVI Fa .A total of 48 ∆VI images were obtained.
Images taken close to the day of the year on which the reference image was taken were selected to reduce the effects of seasonal variation in vegetation.However, the "pseudo-changes" in vegetation caused by differences in the annual data could still be captured by the VI images.The four ∆VI distributions were compared to assess whether such "pseudo-changes" in vegetation caused by changes in the annual data had been captured by the CDR-SR images and whether the Fa-SR images could reduce these effects.
A simple classification extraction was conducted to investigate whether the use of the CDR-SR images, the Fa-SR images and the different VIs (i.e., the NDVI, the GDVI, the GEMI, and the SAVI) affected the mapping of the vegetation changes.Generally, an appropriate threshold was selected to classify the pixels of the ∆VI images as representing a change or no change.Several methods are available to determine these thresholds [34][35][36].Ideally, independent reference data, such as surface data and high-resolution aerial photographs, would be used to verify the accuracy of the classifications [37].However, independent reference data were not available for the study area.The upper and lower 5% were selected as the final thresholds based on multiple experiments that combined the field data, Google Earth and visual maps [38].The upper 5% reflected the positive changes in vegetation, and lower 5% represented the negative changes.The pixels of the 48 ∆VI images were divided into two categories, which reflected either changed or unchanged areas [39].A 3 × 3 median filter was used to reduce the noise of each of the classified images.
The percentage of consistency between the pixels classified as representing the change using the same threshold in the ∆VI CDR images and the corresponding ∆VI Fa images was calculated.A high percentage of consistency would indicate that the method of calibrating the SR had little or no effect on the number of pixels that were classified as reflecting changed vegetation.Moreover, the choice of VI can also affect the number of pixels classified as changed; a high consistency shows that the VI has little effect on the different methods of radiometric correction.

Historical Vegetation Changes
Based on the results of the above analysis, the ∆SAVI Fa images were selected to analyze the historical changes in vegetation in the study area.The annual change of vegetation was not spatially invariant over time [22,40].The changes in vegetation caused by rotation or fallowing of cultivated land or changes in the growth of natural vegetation manifested as steady changes.In the case of wasteland reclamation or serious incursions into previously undisturbed areas, the area of change showed an increasing trend.To investigate these factors, the changed areas (including positive and negative changes) were calculated.
To determine the spatial changes in the study area, six maps derived from the ∆SAVI images reflecting the positive and negative changes were created for all the years to increase the likelihood that all changes in vegetation due to human factors were included.

Comparative Analysis of the VIs
The correlation analysis showed that the four VIs produced using each CDR-SR image displayed a very strong relationship with the corresponding Fa-SR image (R ≥ 0.986).The highest correlation was obtained for the SAVI (R min = 0.991, 1995a), followed by the NDVI (R min = 0.988, 2010a), the GDVI and the GEMI (R min = 0.986, 2010a).The linear regression analysis showed that the largest values of the slope (a) and offset (b) of the four VIs occurred in 1995 compared with the other years.The differences between the pairs of images corresponded to the higher RMSE values in 1995.In general, among the VI images produced from the same CDR-SR images and the Fa-SR images, the lowest values of the RMSE and the strongest correlations were observed for the SAVI images followed by the NDVI images.

Impact of Vegetation Changes
To obtain the actual changes in the surface vegetation rather than spurious effects caused by properties of the images themselves, the areas of vegetation change (including both positive and negative changes) were extracted through a simple classification experiment.The consistency of the classification of the changed vegetation in the ∆VI CDR images was the percentage of changed pixels produced by the CDR-SR images, where the changed pixels were extracted using the CDR-SR images and the Fa-SR images (Figure 3a).In addition, the percentage of changed pixels produced using the Fa-SR images, where the changed pixels were extracted by the CDR-SR images and the Fa-SR images, represented the consistency of changed vegetation in the ∆VI Fa images (Figure 3b).The consistency between ∆VI CDR and ∆VI Fa was poor when the same threshold (5%) was used.The worst consistency was obtained in 1995-2000, followed by 1990-1995, which displayed a significant relationship with the larger slope, offset, and RMSE of the VIs calculated using the two kinds of SR images in 1995.The consistency of the two kinds of SR images in the other periods exceeded 60%.However, the changed pixels produced using the CDR-SR images and the Fa-SR images differed among the different VIs and the different periods.The ∆SAVI images displayed a higher consistency than the other ∆VIs between the CDR-SR images and the Fa-SR images (Figure 3).Fa-SR images, where the changed pixels were extracted by the CDR-SR images and the Fa-SR images, represented the consistency of changed vegetation in the ΔVIFa images (Figure 3b).The consistency between ΔVICDR and ΔVIFa was poor when the same threshold (5%) was used.The worst consistency was obtained in 1995-2000, followed by 1990-1995, which displayed a significant relationship with the larger slope, offset, and RMSE of the VIs calculated using the two kinds of SR images in 1995.The consistency of the two kinds of SR images in the other periods exceeded 60%.However, the changed pixels produced using the CDR-SR images and the Fa-SR images differed among the different VIs and the different periods.The ΔSAVI images displayed a higher consistency than the other ΔVIs between the CDR-SR images and the Fa-SR images (Figure 3).As shown in Figure 3, the worst consistency in the changed pixels contained in the CDR-SR images and Fa-SR images was observed in the periods 1995-2000 and 1990-1995, respectively, and large inconsistencies were observed between the areas of vegetation change obtained using the two types of SR images.Many small changes in the images themselves were captured by the ΔGDVICDR images in 1995-2000 and 1990-1995 (Figures 4 and 5).These changes were monitored for "pseudochanges" in vegetation.These changes were considered to be disruptive to the extraction of vegetation changes and represented the main reason for the poor consistency.Correspondingly, the ΔVIFa images displayed better performance in eliminating the changes in the images themselves and more accurately reflected the changes in vegetation.Most of the areas of the Gobi landscape (Figure 4: Sections A and B; Figure 5: Section F) and bare rock (Figure 4: Section E; Figure 5: Section G) were monitored for "vegetation changes" in the ΔGEMICDR images.The changes in natural vegetation caused by changes in water were also detected by the ΔNDVICDR and ΔGEMIFa images (Figure 4: Section C and D; Figure 5: Section H).The effect of vegetation changes caused by the soil and image brightness factors were eliminated more completely from the ΔSAVI images than from other ΔVI images.The vegetation changes extracted using the ΔSAVI images mainly reflected vegetation changes (either reductions or increases) in the periphery of the protected area.The corresponding ΔSAVI inconsistent maps displayed similarly spatial patterns of inconsistent as the ΔNDVI images; however, the affected area was significantly reduced.As shown in Figure 3, the worst consistency in the changed pixels contained in the CDR-SR images and Fa-SR images was observed in the periods 1995-2000 and 1990-1995, respectively, and large inconsistencies were observed between the areas of vegetation change obtained using the two types of SR images.Many small changes in the images themselves were captured by the ∆GDVI CDR images in 1995-2000 and 1990-1995 (Figures 4 and 5).These changes were monitored for "pseudo-changes" in vegetation.These changes were considered to be disruptive to the extraction of vegetation changes and represented the main reason for the poor consistency.Correspondingly, the ∆VI Fa images displayed better performance in eliminating the changes in the images themselves and more accurately reflected the changes in vegetation.Most of the areas of the Gobi landscape (Figure 4: Sections A and B; Figure 5: Section F) and bare rock (Figure 4: Section E; Figure 5: Section G) were monitored for "vegetation changes" in the ∆GEMI CDR images.The changes in natural vegetation caused by changes in water were also detected by the ∆NDVI CDR and ∆GEMI Fa images (Figure 4: Section C and D; Figure 5: Section H).The effect of vegetation changes caused by the soil and image brightness factors were eliminated more completely from the ∆SAVI images than from other ∆VI images.The vegetation changes extracted using the ∆SAVI images mainly reflected vegetation changes (either reductions or increases) in the periphery of the protected area.The corresponding ∆SAVI inconsistent maps displayed similarly spatial patterns of inconsistent as the ∆NDVI images; however, the affected area was significantly reduced.As showed in Figure 6, obvious changes in the vegetation in the main part of the reserve occurred during 1986-1990; fewer obvious changes occurred in the other periods.These findings indicate that, since 1990, the vegetation in the reserve has generally remained stable.The changed areas were consistently concentrated in the northwest and southwest corners of the reserve, adjacent to the artificial oasis, and were easily affected by human activities.The small plots inside the reserve showed small changes in space and time from 1995-2015.
As is shown in Figure 7, the areas reflecting vegetation changes decreased from 1990 to 1995, increased after 1995, and then decreased again until 2010-2015.The positive changes in vegetation showed the same trend as the overall vegetation changes.The area of positive changes was consistently larger than that of negative changes in the period of 1986-2010, whereas the areas of positive and negative changes were similar during 2010-2015.The area of negative changes in vegetation displayed a slight decrease during 1990-1995, increased after 1995, decreased during As showed in Figure 6, obvious changes in the vegetation in the main part of the reserve occurred during 1986-1990; fewer obvious changes occurred in the other periods.These findings indicate that, since 1990, the vegetation in the reserve has generally remained stable.The changed areas were consistently concentrated in the northwest and southwest corners of the reserve, adjacent to the artificial oasis, and were easily affected by human activities.The small plots inside the reserve showed small changes in space and time from 1995-2015.
As is shown in Figure 7, the areas reflecting vegetation changes decreased from 1990 to 1995, increased after 1995, and then decreased again until 2010-2015.The positive changes in vegetation showed the same trend as the overall vegetation changes.The area of positive changes was consistently larger than that of negative changes in the period of 1986-2010, whereas the areas of positive and

Discussion
The high correlations between the various VI images calculated using the CDR-SR images and the Fa-SR images occurred because the VIs were calculated using a non-linear combination of the red and near-infrared bands of the SR images.The differences in the offset and slope and the RMSE values were caused by the different radiation transmission models employed in producing the CDR-SR images and the Fa-SR images [22,41].The linear regression analysis of the VIs showed that the atmospheric correction method used to generate the CDR-SR images and the FLAASH correction method were not equivalent among the VIs.
Obvious differences in the consistency of the vegetation changes produced using ΔVICDR and ΔVIFa with the same threshold were observed, likely because the CDR-SR images and the Fa-SR images used to calculate the VIs were extracted using different radiation transmission equations [22,42].However, the choice of VIs had a stronger effect on the consistency of the pixels that were classified as changed than the choice of SR images.The ΔGDVI images yielded the worst consistency when the changed vegetation pixels were extracted using the CDR-SR images and the Fa-SR images with the same threshold.Among the ΔVI images, the ΔGDVI images were more sensitive to differences in the images caused by precipitation and different acquisition times; thus, they captured more variable information [14,15].These phenomena were more pronounced in the CDR-SR images than in the Fa-SR images.The ΔGEMI images eliminated some of the differences associated with the images themselves, although the areas of bare rock and Gobi landscape in the study area were also monitored for vegetation changes.However, the changed pixels extracted using the ΔNDVI and ΔSAVI images were found to have little effect on the residual atmospheric correction effect in the images [31,42].Subtle changes caused by water and natural vegetation were also detected by the ΔNDVI images, which reduced the ability of the method used here to distinguish the vegetation changes driven by artificial effects from those produced by natural factors.Since the Donghu Reserve is close to the artificial oasis, and beacuse protective forest management has been implemented there, the vegetation changes were mainly caused by human factors [8][9][10].Thus, the ΔSAVIFa images were used to detect and reflect the changes in vegetation in the Donghu Reserve in this study.
The historical vegetation changes in the Donghu Reserve occurred mainly at the edge of the artificial oasis between 1986 and 2015 (Figure 6).The positive changes in vegetation were consistent with the overall trend in vegetation change (Figure 7).According to Xu et al., Zhang et al. [2,11], and the field investigation, the artificial oasis of the Dang River Plain expanded parallel to the Donghu Reserve.Over the period from 1986 to 1990, the change in the area of the oasis was very small and had little impact on the reserve.The positive changes in vegetation occurred primarily near Xihu

Discussion
The high correlations between the various VI images calculated using the CDR-SR images and the Fa-SR images occurred because the VIs were calculated using a non-linear combination of the red and near-infrared bands of the SR images.The differences in the offset and slope and the RMSE values were caused by the different radiation transmission models employed in producing the CDR-SR images and the Fa-SR images [22,41].The linear regression analysis of the VIs showed that the atmospheric correction method used to generate the CDR-SR images and the FLAASH correction method were not equivalent among the VIs.
Obvious differences in the consistency of the vegetation changes produced using ∆VI CDR and ∆VI Fa with the same threshold were observed, likely because the CDR-SR images and the Fa-SR images used to calculate the VIs were extracted using different radiation transmission equations [22,42].However, the choice of VIs had a stronger effect on the consistency of the pixels that were classified as changed than the choice of SR images.The ∆GDVI images yielded the worst consistency when the changed vegetation pixels were extracted using the CDR-SR images and the Fa-SR images with the same threshold.Among the ∆VI images, the ∆GDVI images were more sensitive to differences in the images caused by precipitation and different acquisition times; thus, they captured more variable information [14,15].These phenomena were more pronounced in the CDR-SR images than in the Fa-SR images.The ∆GEMI images eliminated some of the differences associated with the images themselves, although the areas of bare rock and Gobi landscape in the study area were also monitored for vegetation changes.However, the changed pixels extracted using the ∆NDVI and ∆SAVI images were found to have little effect on the residual atmospheric correction effect in the images [31,42].Subtle changes caused by water and natural vegetation were also detected by the ∆NDVI images, which reduced the ability of the method used here to distinguish the vegetation changes driven by artificial effects from those produced by natural factors.Since the Donghu Reserve is close to the artificial oasis, and beacuse protective forest management has been implemented there, the vegetation changes were mainly caused by human factors [8][9][10].Thus, the ∆SAVI Fa images were used to detect and reflect the changes in vegetation in the Donghu Reserve in this study.
The historical vegetation changes in the Donghu Reserve occurred mainly at the edge of the artificial oasis between 1986 and 2015 (Figure 6).The positive changes in vegetation were consistent with the overall trend in vegetation change (Figure 7).According to Xu et al., Zhang et al. [2,11], and the field investigation, the artificial oasis of the Dang River Plain expanded parallel to the Donghu Reserve.Over the period from 1986 to 1990, the change in the area of the oasis was very small and had little impact on the reserve.The positive changes in vegetation occurred primarily near Xihu nitrate ore, where the vast desert is located, far from the oasis and traffic line and free from human activities; thus, these changes are likely related to natural factors, especially precipitation.The vegetation changes decreased during 1990-1995 because a small increase in the area of the oasis occurred, and the vegetation within the protected area displayed essentially no change during this period.From 1995 to 2010, the artificial oasis near Yingwoliang, Xindiantai, and Dongfeng Farm gradually expanded and encroached on the reserve, resulting in an increase in the area displaying vegetation changes [2].During 2010-2015, because of the small change in the area of the oasis, the area of vegetation change in the study area decreased [11].The slight changes in vegetation in the southeast of the reserve were related to the implementation of enclosure management in the region.Since this area is located far from the artificial oasis, the vegetation changes were less strongly controlled by human factors than the areas near the oasis.vegetation changes in the reserve were mainly related to changes in natural vegetation, and the changes were not obvious.

Conclusions
This study showed that the VI images calculated using the CDR-SR images and the Fa-SR images were strongly correlated (R ≥ 0.986).The differences in the slope and offset were small, but important.The consistency of the vegetation changes detected using the ∆VI CDR and ∆VI Fa images was poor during 1990-1995 and 1995-2000.This poor consistency displayed significant relationships with the larger slope, offset, and RMSE of the VIs calculated using the two types of SR images in 1995.However, the consistency of the vegetation changes extracted using the SAVI images was higher than that obtained with the other VIs, and the ∆VI Fa images displayed better characteristics in terms of the changes in the images themselves compared with the ∆VI CDR images (Figure 3).
The resulting vegetation changes extracted using ∆GDVI CDR and ∆GDVI Fa included many small changes in the images themselves.Most of the areas of the Gobi landscape and bare rock were monitored for vegetation changes in the ∆GEMI CDR images.The changes in natural vegetation caused by changes in water were also detected using the ∆NDVI CDR images and the ∆GEMI Fa images.Thus, the ∆SAVI Fa images were considered to represent the best choice for the extraction of vegetation changes within the study area.
Changes in vegetation have been occurring in the Donghu Reserve since at least 1986.Obvious changes in vegetation in the main part of the reserve were only observed in the period from 1986 to 1990, with fewer changes observed in the other periods.These findings indicate that since 1990, the vegetation in the reserve has generally remained stable.From 1986 to 1990, positive changes in vegetation occurred mainly near Xihu nitrate ore within the protected area.After 1995, the areas of changed vegetation (including positive and negative changes) were mainly distributed near Xindiantai, Dongfeng Farm, and Yingwoliang at the edges of the artificial oasis and along the periphery of the protected area.Positive changes in vegetation dominated the overall trend in the changes in vegetation, with changes decreasing during 1990-1995, increasing until 2005-2010, and then decreasing again during 2010-2015.The trends in vegetation change corresponded strongly with changes in the oasis area of the Dang River Plain, highlighting the increased human pressure on the reserve.In the present study, the analysis of changed vegetation was primarily based on Landsat images.Further applications involving a combination of high temporal and spatial resolution images and hydrological data would be valuable for analyzing changes in vegetation in the Donghu Reserve.

Figure 1 .
Figure 1.The study area is defined by a 2-km buffer that surrounds the Donghu Nature Reserve in Dunhuang (red dotted line).The map on the left is a false-color composite Landsat TM image made up of bands 6, 5, and 4 (i.e., the SWIR1, NIR, and red bands, respectively) for 2015.

Figure 1 .
Figure 1.The study area is defined by a 2-km buffer that surrounds the Donghu Nature Reserve in Dunhuang (red dotted line).The map on the left is a false-color composite Landsat TM image made up of bands 6, 5, and 4 (i.e., the SWIR1, NIR, and red bands, respectively) for 2015.

Figure 2 .
Figure 2. Workflow used in applying the method.

Figure 2 .
Figure 2. Workflow used in applying the method.

Figure 3 .
Figure 3.Comparison of the consistency among the four ΔVIs calculated using the CDR-SR images and the Fa-SR images: (a) the percentage of changed pixels produced using the CDR-SR images, where the changed pixels were extracted using the CDR-SR images and the Fa-SR images; (b) the percentage of changed pixels produced using the Fa-SR images, where the changed pixels were extracted using the CDR-SR images and the Fa-SR images.Horizontal axis labels indicate the pairs of the two years used to create the ΔVI images (centuries are not shown).

Figure 3 .
Figure 3.Comparison of the consistency among the four ∆VIs calculated using the CDR-SR images and the Fa-SR images: (a) the percentage of changed pixels produced using the CDR-SR images, where the changed pixels were extracted using the CDR-SR images and the Fa-SR images; (b) the percentage of changed pixels produced using the Fa-SR images, where the changed pixels were extracted using the CDR-SR images and the Fa-SR images.Horizontal axis labels indicate the pairs of the two years used to create the ∆VI images (centuries are not shown).

Figure 4 .
Figure 4. Maps of inconsistencies in the pixels classified as changed for the period of 1995-2000 in the ΔGDVI, ΔGEMI, ΔNDVI, and ΔSAVI images produced using the CDR-SR images and the Fa-SR images.Annotated regions A to E are referred to in the text.

Figure 5 .
Figure 5. Maps of inconsistencies in the pixels classified as changed for the period of 1990-1995 in the ΔGDVI, ΔGEMI, ΔNDVI, and ΔSAVI images produced using the CDR-SR images and the Fa-SR images.Annotated regions F to H are referred to in the text.

Figure 4 .
Figure 4. Maps of inconsistencies in the pixels classified as changed for the period of 1995-2000 in the ∆GDVI, ∆GEMI, ∆NDVI, and ∆SAVI images produced using the CDR-SR images and the Fa-SR images.Annotated regions A to E are referred to in the text.

Figure 4 .
Figure 4. Maps of inconsistencies in the pixels classified as changed for the period of 1995-2000 in the ΔGDVI, ΔGEMI, ΔNDVI, and ΔSAVI images produced using the CDR-SR images and the Fa-SR images.Annotated regions A to E are referred to in the text.

Figure 5 .
Figure 5. Maps of inconsistencies in the pixels classified as changed for the period of 1990-1995 in the ΔGDVI, ΔGEMI, ΔNDVI, and ΔSAVI images produced using the CDR-SR images and the Fa-SR images.Annotated regions F to H are referred to in the text.

Figure 5 .
Figure 5. Maps of inconsistencies in the pixels classified as changed for the period of 1990-1995 in the ∆GDVI, ∆GEMI, ∆NDVI, and ∆SAVI images produced using the CDR-SR images and the Fa-SR images.Annotated regions F to H are referred to in the text.

Figure 6 .
Figure 6.Changes in vegetation detected in each period from 1986 to 2015.The data of the oasis were described by Zhang et al. [11].Annotated regions (a) to (f) are referred to in the text.

Figure 6 .
Figure 6.Changes in vegetation detected in each period from 1986 to 2015.The data of the oasis were described by Zhang et al. [11].Annotated regions (a) to (f) are referred to in the text.
negative changes were similar during 2010-2015.The area of negative changes in vegetation displayed a slight decrease during 1990-1995, increased after 1995, decreased during 2005-2010, and then increased again.However, the positive changes in vegetation dominated the overall trend of changes in vegetation.Sustainability 2017, 9, 1780 10 of 13 2005-2010, and then increased again.However, the positive changes in vegetation dominated the overall trend of changes in vegetation.

Figure 7 .
Figure 7.The area of vegetation changes during 1986-2015; the horizontal axis tick labels indicate the pairs of the two years (centuries are not shown).

Figure 7 .
Figure 7.The area of vegetation changes during 1986-2015; the horizontal axis tick labels indicate the pairs of the two years (centuries are not shown).

Table 1 .
Landsat acquisition date and sensor type.

Table 1 .
Landsat acquisition date and sensor type.