Spatio-Temporal Dynamics of Land-Use and Land-Cover in the Mu Us Sandy Land, China, Using the Change Vector Analysis Technique

Abstract

The spatial extent of desertified vs. rehabilitated areas in the Mu Us Sandy Land, China, was explored. The area is characterized by complex landscape changes that were caused by different drivers, either natural or anthropogenic, interacting with each other, and resulting in multiple consequences. Two biophysical variables, NDVI, positively correlated with vegetation cover, and albedo, positively correlated with cover of exposed sands, were computed from a time series of merged NOAA-AVHRR and MODIS images (1981 to 2010). Generally, throughout the study period, NDVI increased and albedo decreased. Improved understanding of spatial and temporal dynamics of these environmental processes was achieved by using the Change Vector Analysis (CVA) technique applied to NDVI and albedo data extracted from four sets of consecutive Landsat images, several years apart. Changes were detected for each time step, as well as over the entire period (1978 to 2007). Four categories of land cover were created—vegetation, exposed sands, water bodies and wetlands. The CVA’s direction and magnitude enable detecting and quantifying finer changes compared to separate NDVI or albedo difference/ratio images and result in pixel-based maps of the change. Each of the four categories has a biophysical meaning that was validated in selected hot-spots, employing very high spatial resolution images (e.g., Ikonos). Selection of images, taking into account inter and intra annual variability of rainfall, enables differentiating between short-term conservancies (e.g., drought) and long-term alterations. NDVI and albedo, although comparable to tasseled cap’s brightness and greenness indices, have the advantage of being computed using reflectance values extracted from various Landsat platforms since the early 1970s. It is shown that, over the entire study period, the majority of the Mu Us Sandy Land area remained unchanged. Part of the area (6%), mainly in the east, was under human-induced rehabilitation processes, in terms of increasing vegetation cover. In other areas (5.1%), bare sands were found to expand to the central-north and the southwest of the area.

1. Introduction

The degradation of soil and vegetation in drylands is termed desertification and generally results from combined factors, including climatic variations and human activities. Sandification is a related term that refers to land degradation in sandy areas, involving aeolian erosion, windblown sands, shifting dunes and the encouragement of sands into agriculture and even inhabited areas [1,2]. Certainly, sandification processes are associated with vegetation deterioration, productivity reduction, unstable ecosystems, environmental detriment and economic loss.

Considerable attention has been given to sandification processes in China, since vast areas of sandy deserts are located in the north of the country within arid and semi-arid climatic zones, where the annual precipitation is below 500 mm [1,3]. By the end of 2009, the desertified land area of China was estimated to be 2,623,700 km² (27.33% of the national territory), and the sandified land area to be 1,731,100 km² (69% of the total desertified land area) [4].

Based on literature review, detailed below, the temporal dynamics of the northern China sandy landscapes can be roughly separated into two periods, from 1950 to 1980 and from the 1980s onward. Records show that the rate of development of sandified land in northern China was 1560 km²·yr⁻¹ between the late 1950s and 1975, 2100 km²·yr⁻¹ between 1975 and 1987, and 3600 km²·yr⁻¹ between 1987 and 2000 [5]. During the 1950s, 1960s and the 1970s, the sandification in northern China was the result of interaction between environmental and physical conditions and climatic and anthropogenic factors [6]. Sandy deserts are characterized by sparse vegetation and unconsolidated sands that tend to move under harsh dry and windy weather conditions. At least from the 1970s, this process was intensified due to frequent drought years that were caused by a gradual increase in the annual mean temperature, unstable precipitation, high rates of evaporation, as well as increasing wind speed [7]. Anthropogenic activities were expressed mainly through irrational human economic actions and land use practices, including: (1) overgrazing of about 200% above the carrying capacity with a high ratio of goats to sheep; (2) overcultivation since large-scale grassland was reclaimed into farmland; and (3) unrestricted cutting for the collection of firewood [3,6]. Although it is difficult to effectively separate the impact of anthropogenic desertification from that of environmental desertification due to the insufficient amount of data collected [1], it is believed that the human factors have been the dominant ones in enhancing sandification processes in northern China [5,6,8].

These trends of anthropogenic activities have been reversed since the 1980s. The rate of increasing sandified land has reduced in northern China, mainly due to the great amount of attention paid to this matter by the central government since the end of the 1970s. The state policy was to “plant trees everywhere and make the country green” [9]. This policy was enforced through a series of laws, regulations, resolutions and rules that were then followed by a series of significant strategic measures adopted by the local governments, with close coordination among their agencies, as well as wide public participation and arduous efforts made by local people (e.g., tree-planting campaigns). Actions to combat sandification, undertaken by the government and local residents, including ground and air-seeding of trees, bushes and grasses over large areas, the construction of long windbreaks, shelterbelts and barriers, pastureland enclosures, as well as chemical mulching and hydrologic solutions [2,10]. There is compelling evidence that, during these years, the area of grassland and woodland biomass production enlarged [11–14]. A monitoring program over the sandy deserts of northern China revealed that, between 2000 and 2004, the net area of sandified land decreased annually by 1283 km², and between 2005 and 2009, it decreased annually by 1717 km². In parallel, the average vegetation coverage in the sandified lands rose from 17.63% in 2004 to 18.69% in 2009.

In summary, the above review exhibits that the study area was subjected to different drivers that interact with each other, causing multiple consequences and therefore create complex landscape changes. On one hand, intensive human activities have caused curbing of sandification in some locations with a continuous reduction of the total sandified area. On the other hand, due to changes in climatic conditions (e.g., air temperature, wind speed), observations reveal intensification and expansion in other locations [4,8,15–22]. In addition, it has also been reported that the continued growth in mining and urbanization has had adverse impacts on the landscape [11].

Due to the ability of spaceborne imagery to cover, continuously and repeatedly, vast and remote areas and to evaluate land-cover and long-term changes over the last decades, numerous studies have used remote sensing data to assess the spatial extent of sandification in northern China and its temporal dynamics. This was implemented by various platforms and resolutions along with different image processing methods. Early data, from the 1960s, on land cover were manually interpreted from black and white aerial photographs and supplementary maps [23]. Most studies analyzed Landsat MSS/TM/ETM+ data from selected years on a regional scale [6,8,11–13,15–18,20–22,24–26], supplemented by thematic, administrative and topographic maps. Others used low-resolution imagery data from AVHRR and MODIS [12,21,27]. Using different biophysical indicators (e.g., vegetation indices, fractional cover, albedo, etc.), quantitative and qualitative assessments were applied to classify the severity levels of sandification [13,15,16,18,25,26].

From the above-mentioned studies, one may infer that in the sandy lands of northern China there are two factors that are prevalent and interacted to create different trends of land-cover change. On one side, due to climatic changes, sandification has developed rapidly and is still a serious threat to the fragile ecosystem. However, it is also evident that a considerable portion of the area has undergone rehabilitation processes by means of large cover of grasses, bushes, trees, as well as agricultural crops. Hence, the sandification has been controlled and reversed to some extent. Open questions remain: (1) In a given region, where are the desertified areas and where are the rehabilitated ones? (2) What is the degree of change in each direction? And (3) are these human- or climate-induced changes?

In order to explore these different trends and to consider the knowledge gaps, an appropriate remote sensing based change detection technique has to be employed, taking into account vegetation recovery on one hand, and sandification on the other. In this regard, Change Vector Analysis (CVA) seems to be the most suitable one, since it is a multivariant change detection technique that accepts two or more spectral bands, spectral features or biophysical indicators (e.g., spectral indices, land surface temperature, etc.) associated with the pixel of an image [28,29]. CVA was originally proposed by Malila [28] and further developed by the pioneering studies of Lambin, Strahlers and Ehrlich [30–33]. In CVA, each pixel is represented by a point in a multidimensional space (spectral and temporal) [34]. The methodology (detailed in the Datasets and Methodology section) enables high resolution pixel-by-pixel mapping of the direction and the magnitude of the changes. CVA has several notable advantages over other change detection techniques. Due to its two output products, change magnitude and direction, the CVA technique enables the user to stratify different types of change with respect to land use or land cover dynamics in a specific region. CVA also has the capability to concurrently analyze change in all data layers as opposed to selected bands [34]. Note that most of the studies using this technique have been carried out in forested regions (e.g., [28,35–37]), where the ground features are relatively easily defined (e.g., biomass loss, forest clearing, biomass gain, burnt area, etc.). Far fewer works (e.g., [38–41]) have applied this technique to drylands. In deserts, the biophysical variables and their resulting directional interpretation are different, as will be described later.

Consequently, the primary objective of the current project is to study the land-use/land-cover spatial and temporal dynamics in drylands via an image processing methodology. More specifically, the project aims at utilizing the change vector analysis (CVA) technique using spaceborne imagery over the Mu Us Sandy Land, China. This aim was fulfilled by computing NDVI and albedo from selected Landsat images of 1978, 1990, 1999, 2001, and 2007.

2. Study Area

The Mu Us Sandy Land is one of the most seriously sandified areas in China [3]. It is located in the Ordos Plateau, between l07°20′E and 111°30′E and 37°27′N and 39°22′N, with an area of 38,372 km² (Figure 1). Administratively, the Mu Us Sandy Land lies on the border area of the Inner Mongolia Autonomous Region, Shaanxi Province and the Ningxia Hui Autonomous Region.

The Mu Us Sandy Land has a typical continental semi-arid climate. Annual precipitation ranges from 440 mm in the southeast to 250 mm in the northwest, of which 60%–80% is concentrated in the period from June to August. High inter-annual fluctuations of precipitation characterize the region. Potential evaporation is 2220 mm annually; hence, the aridity index indicates an arid environment. The annual mean temperature is about 6.0 °C–8.5 °C, with monthly mean temperatures of 22 °C in July and −11 °C in January [8]. Over 200 days have a wind speed of more than 5 m·s⁻¹. The topography in the study area gently undulates and slopes from northwest to southeast, with elevations ranging from 1200 to 1600 m above mean sea level. The natural landscape features in the Mu Us Sandy Land include shifting dunes, semi-fixed dunes, fixed dunes, low ridged-land, wetlands, alluvial terraces, lakes, rivers and streams [42]. Most of the land (ca. 77%) is covered by sand with varying thickness [43]. There are several main rivers that cross the southeast of the region and flow into the Yellow River. The ground water table is shallow, 1–3 m below ground, especially in the interdune areas [2]. There are several lakes in the interior of the sandy land, most of which are sodic lakes; others contain chloride, and a few are freshwater lakes [44].

The main vegetation species in sandy grassland (e.g., Stipa bungeana and Poa sphondyloides) cover more than 80% of the sandy land area. Artemisia ordosica and Caragana intermedia are the dominant shrub species; both are drought resistant. The other natural vegetation types include steppe, meadow and shrubs. In addition, there are farmlands, distributed along the river or scattered in the sandy grasslands, and artificial forests and shrubs [8,43]. The grassland in the northwest of the sandy land belongs to Inner Mongolia and is used for grazing only; in the east and south of the sandy land, some grassland areas have been reclaimed into farmland, where both animal husbandry and agriculture are the main economic activities. The sandy land was populated by 774,400 inhabitants in 2000 and 918,000 in 2010. Population density, that was around 10 persons per km² in the 1950s, has continuously increased and reached a level of up to 80 persons per km² in 2002 [45]. According to Zhang et al. [17], between 1960 and 2006, the area of farmlands remained almost at the same level; the areas of woodlands and grasslands increased by 20% and 55%, respectively; the areas of settlements and roads grew 4 times; and the areas of water bodies and unused lands decreased by 57% and 43%, respectively.

3. Datasets and Methodology

3.1. Meteorological Data

Meteorological data for 41 years (1970–2010) were acquired for 5 stations within and in the vicinity of the study area, namely Yulin, Hengshan, Dingbian, Yanchi, and Otog Banner (Figure 1) from the China Meteorological Data Sharing Service System ( http://cdc.cma.gov.cn/). The variables include monthly precipitation, air temperature and wind speed. The annual time series of these three variables is presented in Figure 2, and the departure of precipitation from the normal is illustrated in Figure 3.

3.2. NOAA-AVHRR

Data of the NOAA-AVHRR were obtained from NASA’s Land Long Term Data Record (LTDR) project website ( http://ltdr.nascom.nasa.gov/cgi-bin/ltdr/ltdrPage.cgi). In particular, data were acquired from the AVH09 daily surface reflectance (ρ) product (Version 3). This product is combined from the AVHRR systems onboard NOAA-7, 9, 11, 14, and 16 from 1981 to 1999 in a longitude-latitude geographic projection and at a spatial resolution of 0.05 degrees (ca. 5.55 km² over the study area). For the current project, only data from channel 1 (red) and channel 2 (NIR), and quality assessments for July, August, and September were used. The global image was cropped over the study area. The quality assessment data were used for eliminating clouds and cloud shadows. Daily NDVI was computed according to Tucker [46] as:

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

(1)

The maximum value composite [47] was computed for each pixel for the entire growing season.

Daily albedo (α_AVHRR) values were computed by using Equation (2) after [2]:

$${\alpha }_{AVHRR}=-0.3376{\rho }_{\mathit{red}}^2-0.2707{\rho }_{\mathit{NIR}}^2+0.704{\rho }_{\mathit{red}\square \mathit{NIR}}+0.05256{\rho }_{\mathit{NIR}}+0.0035$$

(2)

Cloud and cloud shadow masks were applied for the daily albedo values that, later, were averaged for the entire growing season.

3.3. MODIS

MODIS data were acquired from the NASA-USGS Land Processes Distributed Active Archive Center (LP DAAC) ( https://lpdaac.usgs.gov/products/). NDVI data were extracted from a 16-Day composite with 500 m of the MODIS Terra and Aqua systems (MOD13A1). The quality assessment information was obtained from the Quality Science Data Set (MODLAND_QA). The strictest criterion was selected, e.g., highlighted pixels with the highest probability to be contaminated with clouds and cloud shadows were excluded from the analysis. Finally, maximum value composites were computed for the three summer months (i.e., July, August and September) for each year. Similarly, albedo values were acquired from the MODIS Terra and Aqua systems (MCD43A3). They are 16-Day composites with 500-m resolution. Yearly data were averaged for the same summer months. For both variables, the time series of the NOAA-AVHRR (1981–1999) was merged with that of MODIS (2000–2010) to produce a continuous series of annual average NDVI and albedo over the entire domain.

3.4. Landsat

Landsat images were selected according to the following criteria: (1) images from the peak of the growing season (either August or September), preferably of the same year; (2) images from years with above-average precipitation; and (3) cloudless images. Due to these strict initial conditions and the fact that the area of interest should be mosaiced by four Landsat images, the number of available images in the archives ( http://glovis.usgs.gov/; http://edcsns17.cr.usgs.gov/NewEarthExplorer/; and ftp://ftp.glcf.umiacs.umd.edu/glcf/Landsat/) was found to be very limited. Consequently, four sets of images, several years apart, were chosen between 1978 and 2007 to fulfill the above-mentioned requirements (Figure 3,

Table 1. Selected Landsat images covering the Mu Us Sandy Land in five time periods.

System	Path/Row *	Path/Row	Path/Row	Path/Row	Denoted in Text
	127/33	127/34	128/33	128/34
Landsat-3 MSS	19-Aug-78	19-Aug-78	20-Aug-78	20-Aug-78	1978
Landsat-5 TM	11-Sep-89	29-Aug-90	15-Sep-88	23-Aug-91	1988
Landsat-7 ETM+	31-Jul-99	01-Sep-00	22-Sep-99	22-Sep-99	1999
Landsat-7 ETM+	6-Aug-02	23-Sep-02	10-Aug-01	10-Aug-01	2001
Landsat-5 TM	12-Aug-07	12-Aug-07	03-Aug-07	03-Aug-07	2007

^*Path and row are according to the Landsat’s World Reference System (WRS)-2. Note that, for the MSS, the path numbers are 137 and 138 rather than 127 and 128, respectively.

). Although all the used images were confined to less than two months in the summer, seasonal variability of the rainfall might also has potential impact on the inconsistency of the data. Bearing this in mind, the raw precipitation dataset was analyzed. No notable change in the seasonal patterns throughout the entire study period was found. A unimodal shape with a significant high peak of rainfall amount is evident in July and August. We conclude, therefore, that during a window of 7 weeks (from 31 July to 23 September), the spectral signal is less affected by seasonal variations and, thus, comparable. In addition to the rainy years, one set of images was used for the same season, but for a drought year (Figure 3, Table 1). Ideally, as in the two edge sets (i.e., 1978 and 2007), all the images were acquired in August of the same year. Both years are characterized by high above-average precipitation. The second set is composed of images from four consecutive years (1988–1991) and denoted hereafter as 1988 for simplicity. The third set is composed of two consecutive years (1999 and 2000) and denoted hereafter as 1999. Lastly, the fourth set is composed of two consecutive years (2001 and 2002) and denoted hereafter as 2001.

Pre-processing of each raw image included one-step radiometric and atmospheric corrections using the dark-object subtraction method [48,49] and the latest radiometric calibration coefficients published [50]. NDVI was calculated as in Equation (1). The albedo values for the MSS were calculated according to Xu et al. [15] (after [51]) and those for the TM/ETM+ according to Liang et al. [52] using the respective equations:

$${\alpha }_{\mathit{MSS}}=0.228{\rho }_{B1}+0.217{\rho }_{B2}+0.183{\rho }_{B3}+{0.373}_{B4}$$

(3)

$${\alpha }_{TM/ETM+}=0.356{\rho }_{B1}-0.130{\rho }_{B3}+0.373{\rho }_{B4}+0.085{\rho }_{B5}-0.0018$$

(4)

where subscript Bn stands for the respective Landsat band number. In each set, four images were mosaiced and cropped with respect to the boundaries of the area of interest (Figure 1). To allow comparison between different sets, NDVI and albedo values were normalized (separately) for the five sets (i.e., 1978, 1988, 1999, 2001 and 2007). This operation is a prerequisite for the change detection algorithm, since the same variables (i.e., NDVI and albedo) are normally distributed with a mean value of 0. Normalization enables the minimization of artifacts (i.e., reducing commission errors) while detecting changes of small magnitudes (i.e., minimizing omission errors) [49,53–55].

3.5. CVA

CVA produces have two outputs—change magnitude and change direction. The magnitude of the change (ΔMG) is computed by the Euclidean distance of a vector between two selected dates, the reference date (t1) and the target date (t2):

$$\left|\Delta MG\right|=\sqrt{{\left(\mathit{Band}{1}_{t2}-\mathit{Band}{1}_{t1}\right)}^2+{\left(\mathit{Band}{2}_{t2}-\mathit{Band}{2}_{t1}\right)}^2}$$

(5)

A threshold, usually in terms of standard deviations, is defined for the magnitude values to distinguish between changed and unchanged pixels. In this study, one standard deviation (STD) was determined as a threshold based on previous experience with change detection analysis [53,56,57]. The direction of the change is computed by the angle (θ) of the change vector from the pixel value at these two dates:

$$\text{tan}\theta =\frac{\mathit{Band}{1}_{t2}-\mathit{Band}{1}_{t1}}{\mathit{Band}{2}_{t2}-\mathit{Band}{2}_{t1}}$$

(6)

This way, the angle depicts the type of change with respect to the mutual change directions within pixels in the two bands. Generally, the number of direction categories is equal to 2ⁿ, where n is the number of spectral bands.

4. Analysis and Results

Long-term analysis of the meteorological variables (Figure 2) reveals that during the 41 years, temperature has significantly increased by about 3 °C, with a sharp step upwards in the late 1990s. Wind speed, on the contrary, was gradually and significantly decreased during the same period. Annual precipitation, however, has remained on the same average level of about 320 mm per year, but with notable fluctuations from year to year (Figure 2). These findings are in agreement with previous studies on the Mu Us Sandy Land [2,7,45].

The merged time series of the NOAA-AVHRR (1981–1999) and MODIS (2000–2010) exhibits continuous series of annual average NDVI and albedo over the entire study area (Figure 4). No significant correlation was found between NDVI and precipitation. However, a significantly high negative correlation (r = −0.63; p << 0.01) was found between NDVI and albedo. NDVI gradually increased with time, and in the same period, the albedo decreased. While increasing NDVI values may resemble the rehabilitation processes of the landscape in terms of fractional vegetation cover and biomass, increasing albedo may be an indicator for sandification processes, since exposed sands have higher reflectance values. These characteristics are the rationale for selecting the NDVI and albedo as the two spectral variables of the CVA for assessing the temporal dynamics of the sandy environment. Consequently, for the current project, NDVI and albedo were used as biophysical indicators (rather than spectral bands). This variable selection results in four directional categories, as schematically described in Figure 5 by the four quadrants of the coordinate system. The first quadrant indicates an increase in both albedo and NDVI, a relationship that represents areas with a significant increase in wetlands. The second quadrant indicates an increase in NDVI and a decrease in albedo, a relationship that represents areas with a significant increase in vegetation. The third quadrant indicates a decrease in both albedo and NDVI, a relationship that represents areas with a significant increase in water bodies, and the fourth quadrant indicates an increase in albedo and a decrease in NDVI, a relationship that represents areas with a significant increase in bare sands.

Based on NDVI and albedo, CVA was applied to: (1) a pair of sequential drought and rainy years (1999 and 2001) as observed in the Landsat images; (2) three pairs of consecutive time periods during the rainy years (i.e., 1978–1988, 1988–2001 and 2001–2007) as observed in the Landsat images; and (3) the entire study period (1978–2007) of the Landsat time series.

In order to assess the short-term climate-induced changes in the regions, CVA was computed on the pair of images that combined the most severe drought years (1999) followed by the most rainy years (2001) (Figure 3). Respectively, rainfall changed from 224 to 428 mm (Figure 2). Figure 6 presents the direction and magnitude results for this pair of images. Since the ground water level is relatively high, intensive rainfall years could add extra water into the biological cycle. Therefore, the most notable change appears as relatively vast patches of natural vegetation in the southeastern part of the study area, and vegetation growth can be recognized in the interdune areas along the longitudinal dunes. These rainy years also flooded the lake in the northwest corner of the image.

Figure 7 illustrates the direction and magnitude results for each of the pairs of consecutive time periods of the rainy years. It can be seen that in the early 1980s, agricultural activities were initiated in the southern part of the study area and along rivers (Figure 7A,B). A small water body, called Hongjiannao Lake, in the north-east of the study area, dried out during this period and exhibits the most pronounced change. In the late 1980s and during the 1990s, sandification was recognized in the center and northern part of the area, while growth of vegetation occurred in the south and especially in the southwest (Figure 7C,D). Magnitude of change ranges from low in the bare sands areas, moderate over the vegetation re-growth areas, and high in the water bodies that were flooded again during this period. Note that the magnitude of the change and the extent of the vegetation patches are less than those observed between 1999 and 2001 due to the post-drought recovery (Figure 6). It can be concluded, therefore, that this change was mainly due to agricultural activity rather than a reaction to rainfall. Between 2001 and 2007, new vegetation was observed in the eastern, northern, and western sections of the area, while more sands appeared in the south in area that was covered before by natural vegetation. In this area, small interdune water bodies dried out and turned to be wetlands. Also, several dams were built along the rivers, storing water for irrigation (Figure 7E,F).

Figure 8 shows the CVA products for the entire study period (1978–2007). In general, it can be seen that vegetation growth and dune stabilization is evident in the eastern Mu Us Sandy Land, while bare sands have expanded to the central-northern and the southwestern parts. These results are consistent with those published in other studies (e.g., [21,26]). Several hot-spot change sites were used for visual interpretation of the results and to confirm the validity of the processes resulted by the CVA variables. This was implemented by zooming-in with very high spatial resolution images acquired by Ikonos, QuickBird or WorldView that are available online on the Google Earth website (Figure 9). Figure 9A (corresponds to Point A in Figure 8A) shows intensive rain-fed agricultural fields in the lowlands. The main crops in the Mu Us Sandy Land are maize and potatoes, sown in May and harvested in October. Figure 9B (corresponds to Point B in Figure 8A) exhibits the Batuwan Reservoir on the Wunding River. It is one of a series of water reservoirs in the region that were built to generate electricity and to supply water for irrigation and industrial purposes [58]. Figure 9C (corresponds to Point C in Figure 8A) demonstrates afforestation activities to stabilize small transverse to barchanoid ridge dunes. Trees, such as Pinus sylvestris var. mongolica, have been planted since the 1980s near the city of Yulin, and bushes, mainly Hedysarum scoparium and H. leave, were seeded from the air [3,8]. Figure 9D (corresponds to Point D in Figure 8A) shows wetlands to the north of Hongjiannao, China’s largest desert lake, which has shrunk by about 10 km² between 1989 and 2007, and its water level is declining by 20 cm annually. The lake is fed by a high level of groundwater and a small river flowing from the north. It is predicted that the lake may vanish in a few decades, due to climate change and human activities, since water has been transported to irrigate the nearby croplands [59,60]. Wetlands are located in the northern part of the lake, near the river inlet, where high albedo values are created by exposed salt and alkali layers, while high NDVI values are caused by vegetation [61]. Figure 9E (corresponds to Point E in Figure 8A) shows Kangjia Bala’er, one of the former water lakes that, since drying up, has become an alkali quarry. Figure 9F (corresponds to Point F in Figure 8A) presents shifting longitudinal dunes that were caused by combined anthropogenic reasons (e.g., grazing pressure) and strong westerly winds that characterize the western part of the Mu Us Sandy Land [21]. These dunes expand into arable lands and urban areas.

Table 2. Change of land-cover area between the studied time periods.

Years	Vegetation (km2)/(%)	Bare Sands (km2)/(%)	Water Bodies (km2)/(%)	Wetlands (km2)/(%)	No Change (km2)/(%)	Total Area (km2)/(%)
1978–1988	1719/4.5	1531/4.0	484/1.3	109/0.3	34,528/90	38,371/100
1988–2001	1397/3.6	1765/4.6	161/0.4	396/1.0	34,651/90.3	38,371/100
2001–2007	2137/5.6	1771/4.6	302/0.8	307/0.8	33,853/88.2	38,371/100
1978–2007	2304/6.0	1950/5.1	461/1.2	310/0.8	33,345/86.9	38,371/100

summarizes the total area changed for each of the studied time periods and for the entire period. It can be seen that the main changes, although small (3%–6% for each category), occur in vegetation-cover exposed sands. Changes of water bodies and wetlands are only of about 1%. The majority of the area (88%–90%) remains unchanged. For the entire period, vegetation cover increased by 6% and bare sands increased by 5.1%. Linking these values with Figure 4, it is suggested that during the study years, vegetation was recovered in some areas, while in others areas, vegetation was removed.

5. Discussion

Various earlier studies have calculated the LULC classification of the Mu Us Sandy Land and elsewhere, using Landsat images in different time steps since the 1970s. These studies presented 3–4 severity categories of desertification (e.g., severe, moderate, slight) or dune stability (e.g., mobile, semi-mobile, semi-fixed, fixed) along with their temporal trends in selected years stemming from the availability of the images [8,13,15,17,25,26]. In order to assess long-term changes, an alternative approach was implemented and described in the current research, namely CVA. CVA has two outputs—direction and magnitude; each of these is presented continuously, pixel by pixel, in the spatial domain. This technique produces a relative change between successive images without the need to struggle with uncertain classification. However, CVA strongly depends on the accurate radiometric calibration of each image, as well as on global (i.e., of all images involved) normalization of the spectral bands or variables.

It is worth mentioning that none of previous studies showed an accuracy assessment since no extensive field surveys were conducted in this vast and remote desert, nor have aerial photographs been continuously available for an objective and independent validation of the classifications [21]. For the same reasons, such a procedure was not possible in the current study. However, as commonly accepted that validation is performed with a higher resolution data than the used one, Ikonos and QuickBird images with 4-m resolution were employed for validation over specific key sites (rather than the entire study domain). Results successfully verified that a reservoir behind a dam is in agreement with the water quadrant of the CVA. The agricultural fields and afforestation areas are in agreement with the vegetation quadrant, the exposed salt and alkali layers are in agreement with the wetland quadrant, and the drying lake is in agreement with the bare soil quadrant.

In this study, the CVA methodology was to build “bottom-up”. First, the temporal dynamics of the two variables, NDVI and albedo, were examined in the Mu Us Sandy Land by utilizing a time series of combined NOAA-AVHRR and MODIS datasets. Upon inferring that the NDVI time series is significantly increasing while that of the albedo is decreasing, the hypothesis about the relationship between these variables and rehabilitation vs. sandification processes was developed.

Change detection studies that rely on irregular satellite images might pose a problem and uncertainty in the results [12,62], mainly due to year-to-year fluctuations of precipitation that can introduce a false-alarm in the long-term interpretation. To overcome this drawback, the long-term LULCC analysis was based explicitly on images of above-average rainy years and from the same phenological stage. Seasonal variability of the rainfall is another factor that might affect the spectral signal, and therefore was carefully examined. Short-term changes, such as a rainy year followed a drought year, reveal changes in different areas than long-term changes. This strengthens the reliability of long-term changes’ results.

In the current study, NDVI and albedo were chosen to represent changes in vegetation and bare sands, respectively. Various previous studies that applied CVA as a change detection algorithm for vegetation and soil used the tasseled cap’s brightness and greenness [35,38,63–67]. It is likely that both biophysical spectral variables would end up with close results. However, the chosen set of variables was found to be more applicative when using a long-term Landsat series. The reason is that there are different units for tasseled-cap transformations derived from different Landsat sensors. Digital numbers (DN)-based transformation values are available for MSS onboard Landsat 1 and 2 and TM onboard Landsat 4 and 5 [68,69]. Needless to say, DN-based transformations on different images might introduce a considerable error due to the different geometries of sun illumination. Moreover, no values are available for Landsat 3 (used in the current study). Reflectance-based transformation values are available for TM onboard Landsat 4 [70] and apparent reflectance-based transformation values for ETM+ onboard Landsat 7 [71]. Obviously, these three types of units make it difficult to process long-term data from different Landsat systems. NDVI and albedo, on the other hand, are calculated with reflectance values.

6. Summary and Conclusions

The spatial extent of desertified vs. rehabilitated areas in the Mu Us Sandy Land, China, and their causes, either natural or anthropogenic, were explored to settle previously-published contradictory findings. In order to fulfill this goal, the CVA technique was implemented using two variables—NDVI and albedo. While changes in NDVI reflect changes in vegetation cover, albedo depicts changes in bare sand cover. Using these variables, four categories of land cover change were computed—vegetation, exposed sands, water bodies and wetlands. The conclusions of this paper are twofold—technical and regional as itemized below.

• Prior to the CVA, high accuracy pre-processing of the raw images, including atmospheric and geometric corrections, as well as normalization of the entire dataset, had to be performed. This was a crucial step without which long-term analyses could not have been performed.

• The main advantage of CVA is the ability to concurrently process the data in the spectral and temporal domains, thus enabling the separation of direction and magnitude of the changes in a continuous manner, rather than in broad qualitative classes, such as slight-, moderate-, or severe land degradation.

• The chosen variables (NDVI and albedo) were selected based on long-term time series analysis of these two variables derived from NOAA-AVHRR and MODIS data. The CVA’s direction and magnitude products detect finer changes compared to separate NDVI or albedo difference/ratio images.

• NDVI and albedo, although comparable to tasseled cap’s brightness and greenness indices, have the advantage of being computed using reflectance values extracted from various Landsat platforms since the early 1970s.

• Each of the four major sectors has a biophysical meaning that was validated in selected hot-spots, employing very high spatial resolution images.

• Careful selection of images, taking into account inter- and intra-annual variability of rainfall, enables differentiating between short-term conservancies (e.g., drought) and long-term alterations.

• Over the entire study period, despite the gradual increase in livestock and consequent growing grazing pressure, some portion of the area, mainly in the east, was under rehabilitation processes, in terms of increasing vegetation cover. Re-growth was a result of increased agricultural activity, planting or air seeding of herbs, grasses, bushes, and trees, and of building lee-shelters, windbreaks and barriers. This is most likely due to the national and local governmental policy to stabilize shifting sands and create pastures by increasing the surface roughness and, thus, decreasing wind speed. In other areas, bare sands expanded to the central-north and the southwest. It is concluded that CVA applied to albedo and NDVI may serve as a useful tool for studying LULC dynamics in arid environment.