Monitoring the Spatiotemporal Dynamics of Aeolian Desertification Using Google Earth Engine

Abstract

Northern China has been long threatened by aeolian desertification. In recent years, all levels of the Chinese government have performed a series of ecological protection and sand control projects. To grasp the implementation effects of these projects and adjust policies in time, it is necessary to understand the process of aeolian desertification quickly and accurately. Remote sensing technologies play an irreplaceable role in aeolian desertification monitoring. In this study, the Zhenglan Banner, which is in the hinterland of the Hunshandake Sandy Land, was considered as the research area. Based on unmanned aerial vehicle (UAV) images, ground survey data, and Landsat images called in Google Earth Engine (GEE), the aeolian desertified land (ADL) in 2000, 2004, 2010, 2015, and 2019 was extracted using spectral mixture analysis. A desertification index (DI) was constructed to evaluate the spatial and temporal dynamics of the ADL in the Zhenglan Banner. Finally, a residual analysis explored the driving forces of aeolian desertification. The results showed that (1) the ADL area in the Zhenglan Banner has been trending downwards over the past 20 years but rebounded from 2004 to 2010; (2) over the past 20 years, the area of slightly, moderately, and severely desertified land has decreased at annual rates of 0.4%, 2.7%, and 3.4%, respectively; (3) human activities had significantly positive and negative impacts on the aeolian desertification trend for 20.0% and 21.0% of the study area, respectively, but not for the rest. This paper explored new techniques for rapid aeolian desertification monitoring and is of great significance for controlling and managing aeolian desertification in this region.

1. Introduction

According to the definition of the United Nations Convention to Combat Desertification, desertification is defined as land degradation in arid, semi-arid, or dry subhumid regions as caused by various factors including climate variability and human activities [1]. Desertification can be divided into aeolian desertification, water erosion, soil salinization, and freezing. In China, most of the desertification belongs to aeolian desertification, which is land degradation caused by the uncoordinated human–land relationship and is marked by sandstorm activities [2]. China is among the countries most seriously affected by aeolian desertification throughout the world. The results of the Fifth Chinese National Desertification Monitoring from 2014 showed that the area of aeolian desertification land (ADL) was 1,721,200 km², which accounts for 17.93% of the national land area [3]. In the face of severe aeolian desertification, all levels of government have made tremendous efforts and a series of ecological projects for sand prevention have been launched such as the Three North Shelterbelt Project, the Beijing–Tianjin Sandstorm Source Control Project, and the Return Grazing to Grassland Project. Many studies [4,5] have shown that the trend of aeolian desertification in several regions have been effectively curbed, the overall ecological situation has been continuously improved, and the ADL area has continued to decrease. Data released by NASA in 2019 show that China has increased the global “greening” over the past 20 years [6].

Remote sensing technologies have become the primary method of desertification monitoring due to the fact of its large-scale data acquisition and continuous observation capabilities as well as the advantages of large and repeated coverage, substantial information, and low cost. The main analysis methods include visual interpretation, supervised classification, and use of single-index or multi-index methods. Visual interpretation is a traditional but widely used method to monitor desertification [7]. It is performed through the establishment of interpretation signs and the comprehensive use of the spectrum, texture, shape, and other characteristics of imagery. However, visual interpretation relies on the interpreter’s understanding of aeolian desertification in the study area, which is time-consuming and has difficulty meeting the quantitative development needs of desertification monitoring. The supervised classification method uses known areas as the basis for image classification and classifies imagery using various classifiers including machine learning algorithms [8,9]. In recent years, more complex deep learning algorithms have been applied in desertification monitoring [10]. However, supervised classification has high requirements for the quantity and quality of training samples, which are prone to problems of poor sample representativeness, high sample subjectivity, and insufficient quantitative evaluation capabilities. Many researchers have chosen a representative indicator to classify ADL, where the most widely used are the fractional vegetation coverage [11], normalized difference vegetation index [12], modified soil adjusted vegetation index [13], and others. However, aeolian desertification is a complex and dynamic process and a single indicator can only reflect changes in one aspect of the desertification, making it difficult to achieve high accuracies. Some researchers have comprehensively used multiple indicators for desertification monitoring including decision trees [14], feature space [15], principal component analysis [16], and the analytic hierarchy process [17]. However, a relatively unified comprehensive indicator system and monitoring method has not yet been formulated.

The causes of aeolian desertification have long been a research focus, and the exploration of the driving forces helps adjust sand prevention and control strategies. Climate change and human activities have been recognized as the two primary driving forces of aeolian desertification. However, accurately and quantitatively analyzing the relative effects of human activities and climate change in aeolian desertification is still difficult. Commonly used methods include qualitative and semi-quantitative assessments [18], regression modeling [19], principal component analysis [20], and residual analysis [21]. Qualitative and semi-quantitative assessments analyze meteorological and social factors based on changes in aeolian desertification and use experienced experts to judge the relative effects. However, this method is subjective and may cause significant differences in the results. The method based on regression modeling analyzes the regression relationship between desertification and climate and social factors. However, the development process of aeolian desertification is complex and nonlinear, making it difficult to simulate using simple regression models. Although the method of principal component analysis has a solid mathematical foundation, it only analyzes the driving forces and does not consider desertification development. In addition, principal component analysis considers the study area without incorporating differences in the continuous space. Residual analysis simulates aeolian desertification conditions without the influence of human factors and uses differences in actual situations to determine the influence of human activities.

The Google Earth Engine (GEE) cloud platform was jointly developed by Google, the United States Geological Survey, and Carnegie Mellon University. The GEE platform stores MODIS data from 2000 to present, Landsat data from 1984 to present, Sentinel data from 2014 to present, and massive amounts of other data such as elevation, meteorology, and land use. The remote sensing data in the GEE platform can be processed in parallel with millions of servers, which significantly reduces the computational time and cost. With its powerful computing capabilities, abundant remote sensing data sources, and convenient processing methods, GEE has advantages in long-term time series and large-scale research and has been widely used to monitor forests [22], agriculture [23], urban areas [24], water [25], etc.

Unmanned aerial vehicles (UAVs) have been developed in recent years as a new type of near-ground remote sensing platform. UAVs have the advantages of low cost and fast acquisition of imagery with high spatial resolutions. They are suitable for acquiring images in areas with complex topographical conditions and that are difficult to access. They effectively compensate for the gap between the regional scale of satellite remote sensing and the sample scale of ground surveys. UAVs can be used for high-precision monitoring of ADL in small areas, which can largely replace traditional manual ground surveys and reduce workloads. To date, UAVs are widely used in vegetation monitoring, habitat monitoring, wildlife protection, etc. [26,27]. However, few studies have applied them to research aeolian desertification [28].

This paper takes the Zhenglan Banner as the research area, which is located within the source of Beijing–Tianjin sandstorms and is in the hinterland of Hunshandake Sandy Land. The Landsat5 and Landsat8 satellite images collected from GEE were used as the primary data resource. Based on the sample points and UAV images obtained in the field, the ADL over the past 20 years was quickly and accurately extracted through spectral unmixing, and a desertification index was established to analyze dynamic changes in the ADL. Finally, the driving forces of aeolian desertification were analyzed.

2. Materials and Methods

The methodological framework of this study is shown in Figure 1, with the key steps of calculation in GEE cloud platform environment. All abbreviations and definitions are shown in

Table A1. Definition of main abbreviations.

Abbreviations	Definition
GEE	Google Earth Engine
UAV	Unmanned aerial vehicle
ADL	Aeolian desertification land
DI	Desertification index
NDL	Non-desertified land
SlDL	Slightly desertified land
MDL	Moderately desertified land
SeDL	Severely desertified land
DSM	Digital surface model
DOM	Digital orthophoto map
SPI	Standardized precipitation index
GLCM	Gray level co-occurrence matrix
MSR	Mobile sand ratio
LSMM	Linear spectral mixture model
RF	Random forest
SMA	Spectral mixture analysis
MNF	Minimum noise fraction
AHP	Analytic hierarchy process

.

2.1. Study Area

The Zhenglan Banner is in the central part of Inner Mongolia and the hinterland of the Hunshandake Sandy Land. It is within the scope of the Beijing–Tianjin sandstorm source control project and is only 180 km away from Beijing. It is approximately 138.7 km long from north to south (41°46’–43°7’N) and 122 km from east to west (114°55’–116°38’ E) with a total area of 10,182 km² (Figure 2). The average annual precipitation of the study area is approximately 365 mm, where the precipitation season is unevenly distributed and concentrated from June to August. The terrain is high in the southwest and low in the northeast with an average altitude of approximately 1300 m. The south part of the study area is mostly low mountains and hills, which is the intersection between the southern edge of Daxinganling and the northern edge of Yanshan Mountains. The north part of the study area is the middle and eastern section of the Hunshandake Sandy Land, which is composed primarily of fixed and semi-fixed sand dunes with a small number of crescent-shaped mobile dunes and low-vegetation coverage.

The vegetation in the Hunshandake Sandy Land is sensitive to climate change and shows the effects of human influence and control, giving it obvious ecological vulnerability [29]. The relatively arid climate and long-term uncontrolled grazing in the study area have made aeolian desertification increasingly serious. The research results on the development status and driving forces of aeolian desertification in the study area will provide a reference for sand prevention projects over the entire Hunshandake Sandy Land and the Beijing–Tianjin Sandstorm Project area.

2.2. Construction of the ADL Classification System

Based on the field survey and previously related research, a land classification system in the Zhenglan Banner was established to include terrain, vegetation, and soil characteristics (

Table 1. Land classification system in the Zhenglan Banner.

ADL Classification	Terrain	Vegetation	Soil
NDL	The terrain is generally flat with open grasslands or lowlands between hills	Psammophytes are companion or occasional species, and the vegetation coverage is greater than 80%	There is no obvious sand-drift activity, the biological crust coverage and silt and clay content are high, and the mobile sand area is less than 15%
SlDL	The topography is not strong and is usually flat inter-dunes or fixed dunes	Psammophytes become the main companion species, and the vegetation coverage decreases to between 50% and 80%	A small amount of sand-drift activity occurs, which is mostly covered by a biological crust or soil layer; the area of mobile sand is between 15% and 40%
MDL	Semi-fixed dunes with small blowouts	Psammophytes become the dominant species, and the vegetation coverage is from 30% to 50%	Sand-drift activity is strong, the biological crust area is greatly reduced, and the mobile sand area is between 40% and 60%
SeDL	The topography is undulating, which includes semi-mobile and mobile sand dunes with large blowouts	There are only psammophytes with a coverage less than 30%	Sand-drift activity is strong, the content of silt and clay particles is extremely low, the mobile sand area exceeds 60%, and the crust surface is small

). The land was divided into non-desertified land (NDL), slightly desertified land (SlDL), moderately desertified land (MDL), and severely desertified land (SeDL).

2.3. Monitoring Period Selection

The five periods of 2000, 2004, 2010, 2015, and 2019 were selected to analyze the aeolian desertification dynamics in the study area. These dates were selected for the following two reasons: First, precipitation is one of the most important factors that affect vegetation conditions and plays an important role in the development of aeolian desertification. Thus, it is only meaningful to analyze changes in desertification under the same precipitation level, otherwise it is difficult to determine land degradation within a certain period. The standardized precipitation index (SPI) from 2000 to 2019 was calculated in R studio by collecting monthly precipitation data for meteorological stations in the study area as shown in Figure 3. We found that the mean SPI values of vegetation growth for 2000, 2004, 2010, 2015, and 2019 were all moderately moist. Second, 2000 was the first year of the Beijing–Tianjin Sandstorm Source Control Project. The five time periods were relatively evenly spaced with a difference of approximately 5 years. Thus, these years represent the overall aeolian desertification conditions over the 20 years of project implementation.

2.4. Data Sets and Preprocessing

This study was conducted at a county level in Zhenglan Banner where medium-resolution remote sensing data was appropriate. The Landsat series of satellites have a long running time, comprehensive image coverage, and rich band information. The surface reflectance products for Landsat5 and Landsat8 are provided in GEE, and the pixel values of each band are the corrected surface reflectance, which can be used directly. In this study, the Landsat5 SR Tier1 and Landsat8 SR Tier1 products were used for the “greenest” composite. First, all images from June to September in the study area were filtered. The clouds and shadows in all images were then masked. Finally, the normalized difference vegetation index (NDVI) of the images was calculated and added to the image as a new band. The “qualityMosaic” function was used to composite all images based on the maximum NDVI, and the annual remote sensing image was obtained by cropping the boundary of the Zhenglan Banner. As this paper focused on the ADL, water bodies, salinized land, and human settlements were extracted and masked in advance to avoid interference from these land-use types.

Precipitation and temperature are important meteorological factors in the aeolian desertification development process. To perform the residual analysis in Section 2.6, this study collected Chinese monthly total precipitation and average temperature images from 2000 to 2017, with a resolution of 250 m. These images were based on the data of more than 1000 meteorological stations and interpolated by AUSPLINE. We selected images of vegetation growth period (May to September) in Zhenglan Banner for further analysis. The MOD13Q1 NDVI data were called in GEE, which was also used for the residual analysis. In addition, statistical data from Zhenglan Banner meteorological station from 2000 to 2019 were collected for analysis.

A field investigation of the Zhenglan Banner was performed along with photographs collected using a UAV in early August of 2019. Global positioning system (GPS) equipment was used to locate various ADL types as classified using Table 1, and a total of 130 sample points were collected. In addition, a same field investigation was performed in the Zhenglan Banner in early September 2010 and 12 sample points were collected. In the follow-up study, these 142 sample points were used to verify the accuracy of ADL classification in Landsat images.

A variety of studies [30,31] have used the DJI Phantom Series UAVs for collection of aerial photographs due to the low cost, simple operation, and mature system. The DJI Phantom4 pro v2 UAV was selected for image acquisition, and 10 representative aerial sample plots with different ADL types were selected (Figure 4). A total of 20 aerial photography sorties were performed. The UAV operated at a height of 80 m with a 75% overlap of the heading and a 65% side overlap. A total of 4382 images were obtained, and the UAV image preprocessing was based on the structure from motion (SfM) method in the pix4dmapper. The digital orthophoto map (DOM) and digital surface model (DSM) were generated, and the pixel size was within 4 cm. Finally, 13 high-quality UAV images were obtained that covered a total area of 1.85 km².

2.5. UAV Image Processing

The UAV images can attain high-precision classification results because of their high spatial resolution. In this study, UAV images were classified as mobile sandy land and non-mobile sandy land. The classification results were then upscaled to obtain the mobile sand ratio (MSR) of the corresponding pixels in the Landsat images, which were used to verify the accuracy of the linear spectral mixture model (LSMM). A total of 11 features were used for the UAV image classification: the original RGB bands, DSM, gray level co-occurrence matrix (GLCM) mean, and 5 vegetation indices (

Table A2. Calculation formulas of UAV vegetation indices.

UAV Indexes	Calculation Formula
Visible band difference vegetation index (VDVI) [48]	$\frac{2green-red-blue}{2green+red+blue}$
Normalized green–red difference index (NGRDI) [49]	$\frac{green-red}{green+red}$
Normalized green–blue difference index (NGBDI) [50]	$NGBDI=\frac{green-blue}{green+blue}$
red–green ratio index (RGRI) [50]	$\frac{red}{green}$
Excess gen index (EXG) [51]	$2green-red-blue$
Color index of vegetation (CIVE) [51]	$0.441red-0.881green+0.385blue+18.78745$

). Before classification, visual interpretation was used to select samples based on the UAV images. A total of 255 samples were obtained, including 131 mobile and 124 non-mobile sand samples. In R studio, the samples were randomly divided into two subsets at a ratio of 1:2 and 75 verification samples, and 180 training samples were obtained. Then, the random forest (RF) model was built based on the “randomForest” package in R studio. This package has two primary parameters to adjust: the number of trees (ntree) and the number of variables (mtry). The grid search method was used to optimize the parameters, showing that the optimal model is obtained when ntree is 600 and mtry is 1. RF can obtain high accuracy classification results when the number of samples is limited, which is suitable for UAV image classification in this study.

A total of 13 UAV images were classified in R studio using the optimal RF model. The method developed by Kattenborn et al. [32] was used to upscale the classification results to match the Landsat pixels. First, the mobile sand pixels were assigned a value of 1 and the non-mobile sand pixels were assigned 0. Second, the results were input into GEE, and the “reduceResolution” command was used to calculate the mean value of it in the corresponding Landsat pixels to obtain the UAV-MSR.

2.6. ADL Extraction Based on Spectral Mixture Analysis

The southern part of the study area is part of the farming pastoral ecotone, where a large area of bare soil appeared from the field survey. Although the vegetation coverage in these areas was relatively low, there was no sand-drift activity, so it could not be classified as ADL. The spectral mixture analysis (SMA) can obtain proportions of various land-use types to give an abundance of each endmember. This method can distinguish mobile sand from other land cover types, and the increase in MSR can represent the enhancement of sand-drift activity and the aggravation of aeolian desertification. In this paper, the range of ADL in the study area was extracted from the endmember abundance of the mobile sand obtained from the LSMM (LSMM-MSR).

In an LSMM, the reflectance of each pixel over the spectral bands is presented as a linear combination of the reflectance of each endmember and its relative abundance, which is defined as

$$\rho \left({\lambda }_i\right)={\displaystyle \sum }_{j=1}^m{F}_j{\rho }_j({\lambda }_i)+\epsilon \left({\lambda }_i\right)$$

(1)

where $j$= 1,2, …, $m$ is the endmember; $i$= 1,2, …, $n$ is the spectral band; $\rho \left({\lambda }_i\right)$ is the reflectance of mixed pixels; ${\rho }_j\left({\lambda }_i\right)$ is the reflectance of endmember $j$ in band $i$; $F_j$ is the abundance of endmember $j$ in the pixel; and $\epsilon \left({\lambda }_i\right)$ is the difference between the actual and modeled reflectance. Field investigations suggested that the three endmembers of mobile sandy land, vegetation, and bare soil can represent the land cover in this area. To find “pure” spectral endmembers, the preprocessed Landsat8 images in 2019 and Landsat5 images in 2010 were output from GEE and were processed using the minimum noise fraction (MNF) rotation in ENVI5.3. Finally, the abundance images of the three endmembers were calculated using the “unmix” command in GEE.

2.7. Dynamic Change Evaluation of ADL

2.7.1. Construction of Desertification Index

To select the most suitable indicators for desertification monitoring and construct the desertification index (DI), 14 indicators (

Table A3. Calculation formula of Landsat indices.

Landsat Indexes	Calculation Formula
Normalized difference vegetation index (NDVI) [33]	$\frac{nir-red}{nir+red}$
Albedo [33]	$\frac{0.356blue+0.130red+0.373nir+0.085swir+0.072swir2-0.018}{1.016}$
Modified soil adjusted vegetation index (MSAVI) [14]	$\frac{2nir+1-\sqrt{{\left(2nir+1\right)}^2-8\left(nir-red\right)}}{2}$
Normalized difference drought index (NDDI) [52]	$\frac{NDVI-NDWI}{NDVI+NDWI}$
Soil moisture monitoring index (SMMI) [53]	$\frac{\sqrt{ni{r}^2+swi{r}^2}}{\sqrt{2}}$
Topsoil grain size index (TGSI) [35]	$\frac{red-blue}{red+blue+green}$
normalized difference soil index (NDSI) [54]	$\frac{swir-nir}{swir+nir}$
Bare soil index (BSI) [46]	$\frac{\left(swir+red\right)-\left(nir+blue\right)}{\left(swir+red\right)+\left(nir+blue\right)}\times 100+100$
$I_{F{e}_2{O}_3}$ [55]	$\frac{\lg \left(red\right)-\lg \left(blue\right)}{\lg \left(red\right)+\lg \left(blue\right)}$

) were initially selected from the aspects of vegetation, soil, humidity, and surface radiation. These indicators were calculated in GEE based on the Landsat image of 2019 (Figure A1). However, selecting too many redundant indicators reduces the accuracy and efficiency of aeolian desertification monitoring. A total of 130 samples obtained from a field investigation in 2019 were used for the gain ratio (Figure A2) and RF (Figure A3) feature selection, and a correlation analysis (Figure A4) was performed for all indicators. Finally, 4 indicators (albedo, NDVI, wetness, and TGSI) with high importance and low correlation were selected to construct the DI.

The albedo is one of the most frequently used indicators in ADL monitoring [33,34] and represents the ratio of the solar radiation flux that is reflected by the surface to the incident solar radiation flux. The topsoil grain size index (TGSI) was proposed by Xiao et al. [35], which reflects the mechanical composition of the topsoil and is positively correlated with the fine sand content. Land surface humidity also has a strong correlation with the degree of aeolian desertification. The wetness of the tasseled cap transform can reflect the humidity status of the soil and vegetation, which has been often used to monitor ecological environments [36]. The NDVI is sensitive to vegetation changes and can reflect the growth status of vegetation on the ADL.

The analytic hierarchy process (AHP) has been widely used in the study of aeolian desertification [37,38] and can reveal the relative importance of different desertification indicators based on the study area. The weight of each variable in the AHP is determined from the matrix of pairwise comparisons, where the importance of indices is ranked from 1 to 7, otherwise they are assigned 1/2 to 1/7 (

Table 2. Relative importance of two indicators.

1/7	1/6	1/5	1/4	1/3	1/2	1	2	3	4	5	6	7
Extremely low		Strong		Moderate		Equal	Moderate		Strong		Extreme important

). Based on many relevant experts and previous research results [17], the final comparison matrix is as shown in

Table 3. Index comparison matrix.

	Albedo	NDVI	TGSI	Wetness	Weights
Albedo	1	2	3	5	0.482
NDVI	1/2	1	2	3	0.272
TGSI	1/3	1/2	1	2	0.157
Wetness	1/5	1/3	1/2	1	0.088

.

As the dimensions of each index are not unified, they should be standardized before any analysis. The standardization methods for the NDVI sand wetness are shown in Equation (2), those for the TGSI and albedo are shown in Equation (3), and the DI is shown in Equation (4):

$$P_i=\frac{max{X}_i-X_i}{max{X}_i-min{X}_i}$$

(2)

$$P_i=\frac{X_i-min{X}_i}{max{X}_i-min{X}_i}$$

(3)

$$DI={\displaystyle \sum }_{i=1}^n{W}_i\times P_i$$

(4)

where $X_i$ represents the $i$th index, $P_i$ is the standardization value of $X_i$, $W_i$ is the weight of the $i$th indicator, and $DI$ is a dimensionless value. A larger $DI$ gives a higher degree of aeolian desertification and a worse ecological status.

2.7.2. Classification and Dynamic Change in ADL

The ADL extracted in Section 2.4 was classified into SlDL, MDL, and SeDL based on the DI. The accuracy of the classification results was verified using 142 samples obtained from the field investigation, and the area of each ADL type was counted in GEE. In addition, the annual change rate for each type of ADL was calculated as

$$\mathsf{\Delta }U=\frac{U_{\left(b,i\right)}-U_{\left(a,i\right)}}{U_{\left(a,i\right)}}\times \frac{1}{T}\times 100\%$$

(5)

where $\mathsf{\Delta }U$ is the change rate of a land type, $U_{\left(a,i\right)}$ is the area of type $i$ at the beginning of a certain period, $U_{\left(b,i\right)}$ is the area of type $i$ at the end of a certain period, and $T$ is the duration of the study period.

To further analyze the spatiotemporal dynamic change of the ADL in the study area over the past 20 years, the change patterns were classified into significantly reserved, reserved, stable, developed, and seriously developed. Reserved is the desertification degree in a certain pixel decreased by a level. Significantly reserved is the desertification degree in a certain pixel decreased by two or more levels (e.g., a change from SeDL to SlDL). Stable is the desertification degree in a certain pixel with no change. Developed is the desertification degree in a certain pixel increased by one level. Seriously developed is a certain pixel with a cross-level increase in the degree of desertification.

2.8. Residual Analysis

Residual analysis judges the effects of human activities on aeolian desertification. The impact of human activities is reflected primarily on vegetation. Based on previous studies [21,39], this paper selected the NDVI to characterize the effects of human activities on aeolian desertification. The mean of the NDVI in the growing season from 2000 to 2017 was calculated in GEE. The regression model for the meteorological factors (average temperature and total precipitation during the growing season) and the NDVI was established using the random forest algorithm in GEE to obtain the NDVI prediction ($NDV{I}_{predict}$). The $NDV{I}_{predict}$ represents the aeolian desertification as affected only by climate. The calculation method of $Residua{l}_{NDVI}$ is shown in Equation (6), which represents the degree of influence of human activities on aeolian desertification.

$$Residua{l}_{NDVI}=NDVI-NDV{I}_{predict}$$

(6)

The Sen slope of the $Residua{l}_{NDVI}$ from 2000 to 2017 was then calculated, and the Mann–Kendall (M–K) significance test was performed. This paper divides human activities into significant positive, insignificant, and significant negative effects based on the impact on the ecological environment. When the Sen slope is >0 and p < 0.05 in the M–K test, it was considered a significant positive effect, which indicates that human activities positively influence ecological conditions and may promote the reversal of aeolian desertification. When p > 0.05, it was considered an insignificant effect, indicating that human activities have no meaningful impact on aeolian desertification. When the Sen slope is <0 and p < 0.05, it was considered as a significant negative effect, indicating that human activities negatively influence ecological conditions, and unreasonable land use may promote aeolian desertification.

3. Results

3.1. UAV-Based Mapping of MSR

Based on the input features and the RF model after parameter adjustments, the overall classification accuracy of 13 UAV images reached 97.6% and the kappa coefficient reached 0.95. The high-precision classification results helped to obtain accurate UAV-MSR samples (Figure 5). After manually deleting pixels with abnormal values and incomplete coverage, a total of 2428 pixels were used to determine the UAV-MSR samples, which was then used to verify the accuracy of the LSMM. Compared with selecting samples from the field and visually determining the MSR in previous aeolian desertification studies [40,41], UAVs provide a larger number of more objective and accurate samples with faster acquisition speeds. This effectively compensates for the scale gap between the samples and remote sensing image pixels.

3.2. Extraction of the ADL

3.2.1. Accuracy Analysis

Based on the UAV-MSR sample pixels obtained in Section 3.1, we obtained the corresponding LSMM-MSR pixels for the Landsat images in GEE. Then, we analyzed the LSMM-MSR based on the UAV-MSR to verify the accuracy of the LSMM. The regression analysis results show that the two data sets had a good linear relationship with an R² value of 0.7648 and RMSE of 0.1383 (Figure 6). This indicates that the LSMM accurately obtained the MSR.

The LSMM-MSR was used to extract the ADL. The experiments showed that the ADL of the study area could be accurately extracted when the LSMM-MSR was 10–20%. To further determine the optimal threshold, 142 samples obtained from field surveys in 2019 and 2010 were used to verify the classification accuracy as shown in

Table 4. Extraction accuracy of the ADL.

LSMM-MSLA	10%	12%	14%	16%	18%	20%
Correctly classified samples	108	115	120	127	129	121
Accuracy	76.1%	81.0%	84.5%	89.4%	90.8%	85.2%

. The highest classification accuracy was obtained when the LSMM-MSR was greater than 18%. Based on this threshold, the preprocessed images were used to extract the ADL of the study area in GEE.

3.2.2. Extent of the ADL

The ADL extraction results of the study area from 2000, 2004, 2010, 2015, and 2019 are shown in Figure 7. ADL was mainly distributed in the north of Zhenglan Banner and concentrated in the eastern and western parts. Over the past 20 years, the ADL area in the east decreased significantly, and now it is mainly concentrated in the northwest. Statistical results of the ADL, shown in Figure 8, indicate a decreasing trend in the study area over the past 20 years. In 2000, the ADL reached 5161.1 km², which accounts for 50.5% of the study area. From 2000 to 2004, the ADL rapidly decreased to 3097.5 km², and its proportion decreased to 30.3%. From 2004 to 2010, the desertification showed a rebounding trend, and the ADL grew to 4462.9 km², which accounts for 43.7%. After 2010, desertification continued to reverse. In 2015, the ADL decreased to 3736.6 km², which accounts for 36.6%, and in 2019, it decreased to 3179.7 km², which accounts for 31.1%. Although there were some fluctuations, the desertification trend in Zhenglan Banner was noticeably reversed.

3.3. Aeolian Desertification Dynamics from 2000 to 2019

Based on the results in Section 3.2.2, the NDL was masked and the DI was used to evaluate the ADL dynamics in the study area. The DI results are shown in Figure 9.

3.3.1. Classification and Time Variation Characteristics of ADL

According to experiments, ADL could be classified into 3 types by DI: SlDL (DI ≤ 0.35), MDL (0.35 < DI ≤ 0.41), and SeDL (0.41 < DI). The confusion matrix of the four land types (

Table 5. Classification accuracy of ADL.

	NDL	SlDL	MDL	SeDL
NDL	39	9	0	0
SlDL	4	27	3	2
MDL	0	4	21	4
SeDL	0	1	3	25
Overall accuracy	78.9%

) was obtained, showing an overall accuracy of 78.9%. The classification results over the considered years are shown in Figure 10.

Figure 11 shows the proportions of various land types. The SeDL area decreased from 1410.2 km² in 2000 to 497.8 km² in 2019, with a decrease of 9.0% of the total area. From 2010 to 2015, the SeDL area decreased by 5.4%, and the restoration of aeolian desertification was significant. The MDL area exhibited a decreasing trend; however, this area surged by 8.1% from 2004 to 2010. The SlDL area decreased from 1977.1 km² in 2000 to 1810 km² in 2019; it has changed little but continued to increase from 2004 to 2015.

The annual change rates (Figure 12) of SlDL, MDL, and SeDL from 2000 to 2019 were –0.4%, –2.7%, and –3.4%. The annual change rates of SeDL for the four periods were –5.5%, 0.2%, –10.0%, and –2.6%. It was the ADL type with the fastest decline over the 20 years. The annual MDL change rate fluctuated greatly; it decreased annually by 12.4% from 2000 to 2004 and increased annually by 15.3% from 2004 to 2010. The annual rates of SlDL were –11.1%, 8.0%, 6.6%, and –4.2% for the four periods. The three types of ADL showed a downward trend over the past 20 years, especially the areas of SeDL and MDL dropped rapidly. This indicates that the aeolian desertification of Zhenglan Banner reversed and the grassland vegetation status achieved marked improvement.

3.3.2. Spatial Variation Characteristics of ADL

The spatial dynamics of the ADL in the study area are shown in Figure 13. The statistics in Figure 14 show that the area of stable during the two periods is the largest, which suggests that the change in aeolian desertification was relatively smooth and there were no particularly drastic changes. Between 2000 and 2004, there was a strong reverse trend in aeolian desertification, and the areas of reversed and significantly reversed land were larger than in other years. The areas of developed and seriously developed land were smaller than in other years and were concentrated mostly in the eastern region. From 2004 to 2010, both the developed and seriously developed land areas were significantly larger than those in other years, which took up the majority of the image as illustrated in Figure 13. The reversed area was distributed mostly in the original SeDL in the northwest and east. From 2010 to 2015, the reversed area accounted for 20.9% of the total, and the significantly reversed area accounted for 6.5%. Few developed areas were distributed in the middle of the study region. From 2015 to 2019, the stable area was larger than in other years, and aeolian desertification still showed a reversed trend. However, there were some developed areas in the northwest. In general, the total proportion of reversed and significantly reversed land from 2000 to 2019 reached 37%. The reverse of the original SeDL area in the east and northwest was the most significant. In addition, the total proportion of developed and seriously developed land was only 5.4% and was scattered at the eastern and western edges.

3.4. Impact of Human Activities on Aeolian Desertification

Figure 15 shows that areas with significant positive effects of human activities were distributed primarily in the middle eastern and western parts of the study area. These areas have the most serious aeolian desertification and are the significant reversal areas, showing that various desertification control projects have played important roles in ecological restoration. Some areas in the south also showed significant positive effects from human activities. The visual interpretation indicates that these areas were newly increased cultivated land over the past 20 years, and the vegetation status was improved from crop planting, which also explains the rationality of the residual analysis. Human activities in the eastern and northwestern parts of the study area have significant negative effects that may be related to overgrazing. Thus, increased human activities in these areas played a role in promoting desertification. The significant negative effects in the south are mostly related to the decreased cultivated land area.

The statistical results showed that the area with insignificant effects accounted for 59.1% of the total, indicating that human activities had no serious impact on aeolian desertification on more than half of the land. The proportion of the area with significant positive effects was 20.0%, and the proportion of significant negative effects was 21.0%. This indicates that the area with positive and negative effects of human activities in the study area were nearly the same over the past 20 years. With the development of various ecological projects, the original SeDL had significantly improved, but there are still large areas that suffer from the unreasonable human activities.

4. Discussion

In this study, a ground survey, UAV images, and satellite images were used to analyze the dynamic of aeolian desertification. Data from three different scales were effectively integrated. In other studies of vegetation coverage [32], the satellite–aviation–ground integration method was often used, but it was rarely used in desertification monitoring. In addition, MSR is more effective than vegetation coverage in extracting ADL. Some existing studies in the Zhenglan Banner used the vegetation coverage [42,43] to classify the ADL, which gave incorrect classifications for most of the bare soil in the south as SlDL or even MDL. Moreover, a multi-index method can more reasonably reflect the degree of ADL. As desertification is a complex process, a single index can only show one aspect of change [17]. The methods used in this paper can provide reference for rapid and accurate monitoring of aeolian desertification in other areas.

The results showed that the area of ADL and the degree of desertification in the study area decreased from 2000 to 2019. This is similar to many related studies in recent years. Yu et al. [44] studied the vegetation change trend in the Beijing–Tianjin Sandstorm Source Control Project area, and the results showed that 83.3% of the vegetation belonged to restoration state from 2005 to 2015. Moreover, the monitoring results of the desertification process in Hunshandake Sandy Land by Gou et al. [45] were basically consistent with those in this paper. As the reversal trend of aeolian desertification is obvious in this area, it is important for us to judge the most important driving factor behind.

The results of the residual analysis show that the aeolian desertification was not significantly affected by human activities for 59.1% of the land area in the study region. This means that the dynamic of aeolian desertification in most areas was affected by climate. Statistical data (Figure 16, Figure 17 and Figure 18) from Zhenglan Banner meteorological station in the past 20 years show that the annual precipitation and average temperature present an upward trend. The warm and humid climate was beneficial to the reversal of aeolian desertification [46]. Although the average annual wind speed increased significantly from 2005 to 2011, it showed a downward trend over the entire 20 years, which reduces wind erosion. Thus, we believe that climate is the most important factor in the change of aeolian desertification in the study area, which is consistent with recent relevant studies. For example, Yang et al. [47] considered the meteorological factor in Hunshandake Sandy Land was the key factor that affects aeolian desertification, while human activities may be the most important factor in other sandy lands throughout China. The increase in ADL area in 2010 may be related to the arid climate in 2009. The total precipitation in 2009 was the lowest in 20 years. At the same time, the annual average wind speed also reached the maximum in 2009 and 2010.

Human activities also had an important impact on aeolian desertification in the study area. The proportion of the total land area with significant positive effects from human activities was 20.0%, most of which was in the SeDL. This indicates that since 2000, sand control measures, such as aerial seeding and grazing prohibition, have played an important role in promoting the reversal of SeDL. However, 21.0% of the area was still affected by the significant negative effects of human activities, which indicates that there are still large areas of damaging land-use practices in the study area. It is still necessary to reduce unreasonable human activities in the future and adjust resource utilization to maintain the restoration of ecological environment and further reduce the ADL area.

5. Conclusions

This paper took the Zhenglan Banner as the research area and explored the processes and driving forces of aeolian desertification from 2000 to 2019. The main conclusions are as follows:

(1) UAVs can quickly provide the MSR. In this paper, the UAV images were divided into mobile sandy land and non-mobile sandy land with an accuracy that reached 97.6%. The classification results were used for resampling in GEE, and a total of 2428 UAV-MSR samples corresponding to the Landsat pixels were obtained to verify the LSMM accuracy.

(2) The ADL area in the study area decreased over the past 20 years. The ADL area decreased from 5161.1 km² in 2000 to 3179.7 km² in 2019 but increased by 1365.4 km² from 2004 to 2010.

(3) Over the past 20 years, the aeolian desertification of the study area showed a reversing trend. The area of SlDL, MDL, and SeDL decreased at annual rates of 0.4%, 2.7%, and 3.4%. In terms of spatial dynamic changes, the total proportion of reversed and significantly reversed lands from 2000 to 2019 reached 37%, and the total proportion of developed and seriously developed lands was only 5.4%.

(4) Climate change and human activities have jointly promoted the reversal of aeolian desertification in the study area, where climate factors played a major role. The results of residual analyses showed that aeolian desertification was positively affected by human activities over 20.0% of the land area, and 21.0% was affected by the negative effects of human activities.

In order to further control aeolian desertification in the study area, it is still necessary to reduce unreasonable human activities and increase project investment in SeDL areas.