Land Cover Mapping and Ecological Risk Assessment in the Context of Recent Ecological Migration

Abstract

In order to protect the ecological environment and solve the poverty problem in the western region, China has established an ecological migration (EM) policy. This policy aims to relocate populations from poverty-stricken areas with fragile ecological environments, which inevitably leads to changes in land cover and the ecological environment. The objective of this study was to identify the effects of EM in a typical region (Wuwei), including changes in the land cover and ecological risk (ER). A land cover change monitoring method was implemented for the 2010–2019 period for six land cover classes using random forest, which is an effective supervised machine learning method. The land cover change patterns were analyzed by determining the area changes of the six classes and applying a land use transition matrix, and a landscape ecological risk model based on landscape disturbance and fragility was used. Our results demonstrate that the increase and decrease in the area of cultivated land, unused land, and construction land can be divided into two stages (2010–2015 and 2015–2019). The area of water and perennial snow doubled during the study periods. The major land cover transitions were between unused land and construction land and between unused land and crop land. In addition, the ER value for the Qilian Mountain National Nature Reserve decreased because of the implementation of EM in the study area, indicating that the ecological environment was effectively improved. The results demonstrate the advantage of the proposed approach in understanding the impact of EM on regional land cover changes and the ecological environment so as to provide guidance for follow-up planning and development.

1. Introduction

Poverty eradication and ecological environment protection are always important objectives of sustainable development [1]. However, the poverty rates are high and the ecological environment is fragile in regions of refugee habitation [2]. In order to effectively solve these two problems, the governments of some countries have proposed a series of remedial measures in which ecological migration (EM) is a key component [3,4,5,6]. Experts of the United Nations Population Development Fund have reported that there are tens of millions of environmental refugees in the world who need to migrate [7]. This project aims to migrate the population from an ecologically vulnerable area to improve the service function and value of the ecosystem and to develop the economy [8]. EM plays an important role in achieving the goal of eliminating absolute poverty [9,10]; however, EM has caused dramatic changes in land cover and the ecological environment. Monitoring the annual change in land cover and analyzing the spatial–temporal patterns of ecological risk (ER) are essential to support regional sustainable development.

EM is a typical human activity and has generated considerable recent research interest [11,12]. Remote sensing techniques are effective instruments to monitor the impact of human activities on land cover changes. Previous studies have focused on changes caused by agricultural land abandonment and retirement [13,14], ecological restoration programs [15,16], greenhouses [17], urban sprawl [18], and so on. In contrast to these activities, the changes resulting from EM include two aspects: building the settlement area and restoring the landform of the emigration area. Therefore, the land cover changes are more complicated. Among the various methods used to detect land cover changes, the approach of first monitoring and then classifying has been widely used [19,20]. This method has proven to be efficient, but it requires the high performance of hardware to process a large number of remote sensing images for multi-temporal land cover change monitoring (LCCM). Fortunately, the explosive growth of Earth observation data and the rapid advancement of cloud computing technology have made large data calculations possible [21]. Google Earth Engine (GEE) is the most representative platform, which was launched by Google in December 2010. This tool provides access to petabytes of publicly available remote sensing imagery, ready-to-use products with an explorer web app, high-speed parallel processing, and machine learning algorithms using the computational infrastructure of Google. Therefore, GEE has led to great achievements in the field of remote sensing and has attracted widespread interest for land cover classification, LCCM, and the extraction of specific objects [22,23,24,25,26].

Ecological changes are a primary manifestation of human utilization of the natural environment. A clear understanding of changes in the ecological environment resulting from human activities will help in formulating and managing sustainable development policies [27,28]. The landscape ecological risk index (ERI) is an important indicator of ecological security. This indicator primarily depends on the coupling of landscape patterns and ecological processes to achieve the integrated characterization of multi-source risks from natural or human activities [29,30]. Researchers have previously explored the impacts of road network expansion [18], geomorphological regionalization [31], and topography [32] on ER. Existing research has focused on the motivation of EM [33], the relationship between EM and the environment [12], and the satisfaction of immigrants [34]. There is currently a lack of research that analyzes the ER changes resulting from the EM project.

In summary, our aim was to describe the spatiotemporal land cover change types and the condition of the ecological environment with the help of a multi-temporal LCCM method and the ERI. Therefore, this study detected land cover changes in Wuwei region from 2010 to 2019 and then evaluated the effects of EM. The specific research objectives were to (1) generate land cover maps for the period from 2010 to 2019, (2) map spatiotemporal changes in the study area during the same period, and (3) quantify the ER from 2010 to 2019.

2. Study Area and Dataset

2.1. Study Area

The study area, Wuwei, lies in the middle of Gansu Province, China (as shown in Figure 1). This prefecture-level city covers an area of 32,347 km². It is at the intersection of the Loess Plateau, the Qinghai Tibet Plateau, and the Mongolian New Plateau [35]. The hypsography is complex and can be divided into three zones. The southern zone is the Qilian Mountains, with an altitude of 2100–4800 m. The climate is cool with abundant precipitation, which is conducive to the development of forestry and animal husbandry. The central zone is a plain oasis area with an altitude of 1450–2100 m above sea level. With flat terrain and fertile land, it is an important production base of grain, oil, fruits, and vegetables for the whole province and even the country. The northern zone is a desert area with less precipitation.

The south of Wuwei is part of the Qilian Mountain National Nature Reserve, which is an important ecological security barrier in the west of China; however, the ecological environment of this area is fragile [32], and the proportion of the population in poverty was high. Therefore, the Wuwei Municipal Party Committee proposed and implemented a large-scale EM project in 2011 to gradually move 72,000 poor farmers and herdsmen living in the high and deep mountain area to a low altitude area. By the end of 2019, the initial effect had been attained, which can be monitored by remote sensing. Therefore, the selection of this study area is of representative significance.

2.2. Dataset

As previously mentioned, the experiment on land cover classification was carried out using GEE. This platform provides online access to archived Landsat data as a collection of the United States Geological Survey (USGS) [36]. We chose Landsat data in their atmospherically corrected surface reflection forms, including Landsat 5 TM from 2010 to 2012, Landsat 7 ETM+ from 2010 to 2019, and Landsat 8 OLI/TIRS from 2013 to 2019, on GEE (in Supplementary Materials code 1). We used the seasonal classes of the Northern Hemisphere (Spring: March to May; Summer: June to August; Autumn: September to November; Winter: December to February). The Landsat images available for each season in 2010–2019 are shown in Figure 2. Image preprocessing included the following steps. Firstly, the Landsat images of the study area in 2010–2019 were selected. Secondly, the CFmask algorithm [37] was applied to mask clouds and create cloud-shadow-free images. Thirdly, to compensate for data loss, we used the following two measures. First, a time series of Landsat image composites in each season was reduced by calculating the median of each pixel across the stack of all matching bands. Second, except for 2012, Landsat 5 TM and Landsat 8 OLI/TIRS images were obtained as the baseline data, and Landsat 7 ETM+ images were used to compensate for data loss. Landsat 5 TM and Landsat 8 OLI/TIRS images were unavailable for 2012, so only Landsat 7 ETM+ images were used. SLC-off images (when the Scan Line Corrector failed and Landsat 7 ETM+ products had data gaps) were processed by calculating the median value in each quarter. In addition, because there were different land cover types at different altitudes, Shuttle Radar Topology Mission (SRTM) data were used [38].

3. Methods

The workflow of the method had four steps: classification of images from the basic year, detection of potential area changes in the comparison years, mapping of annual classification results, and impact assessment of EM (Figure 3). We chose 2019 as the basic year and 2010–2018 as the comparison years. Random forest (RF) was selected as the classifier in this study [39].

3.1. Reference Dataset for Training and Validation Samples

Samples have a crucial impact on classification results [40]. Therefore, selecting representative samples for model training and accuracy evaluation is an important prerequisite to ensure the accuracy of the final map. A training dataset was collected to train the model, while a validation dataset was constructed to test the accuracy of the classification results. Samples were acquired with the help of high-definition resolution images in Google Earth through visual interpretation. As a result, dozens of polygons were selected for each land cover class for each year. Polygons for each year were randomly divided into the training and validation datasets in a proportion of 3:1. In order to prevent imbalanced datasets from affecting the accuracy of minority instances, samples of each class were selected equitably. Then, 1500 points of each class in the training dataset were randomly selected, and 500 points of each class within the validation dataset were randomly collected for accuracy verification. The land cover types were unused land, water and perennial snow area, cropland, forest land, grassland, and construction land.

3.2. The Classification of Basic Year Images

As a result, because the land cover map of 2019 was used as the baseline for further updates of the land cover maps of comparison years, a precise classification result of 2019 was crucial for subsequent stages. In this step, we comprehensively considered various features in order to ensure that the information contained in remote sensing images could be fully detected, including quarterly and yearly spectral features and quarterly textural and terrain features. In addition to six spectral bands of Landsat data (blue, green, red, near infrared, shortwave infrared 1, and shortwave infrared 2), we chose the Normalized Difference Vegetation Index (NDVI) [41], the Enhancement Vegetable Index (EVI) [42], the Green Normalized Difference Vegetation Index (GNDVI) [43], and the Soil-Adjusted Vegetation Index (SAVI) [44] to distinguish vegetation from other instances. We selected the Normalized Difference Build Index (NDBI) [45] to detect construction land and the Modified Normalized Difference Water Index (MNDWI) [46] to distinguish water. The Gray Level Co-occurrence Matrix (GLCM) [47] of the near-infrared band was calculated with a moving window size of 5. Textural features were used to distinguish different land cover types and avoid commission of the same type. Six textural features with low correlations were selected, including angular second moment (ASM), contrast (CON), correlation (COR), variance (VAR), inverse difference moment (IDM), and sum entropy (SENT). Because land cover types are greatly affected by topography, we calculated the terrain elevation, the slope, the aspect, and the hill shade from Shuttle Radar Topography Mission (SRTM) data. The features used in this study are shown in

Table 1. Features in classification.

Spectral Features	Season (*41)	B/G/R/NIR/SWIR1/SWIR2/NDVI/EVI/GNDVI/SAVI/NDBI/MNDWI 2
	Year	NDVImax/NDVImin/EVImax/EVImin/GNDVImax/GNDVImin/SAVImax/ SAVImin/NDBImax/NDBImin/MNDWImax/MNDWImin
Textural Features	Season (*41)	ASM/CON/COR/VAR/IDM/SENT 3
Topographical Features	Year	Elevation/Slope/Aspect/Hill shade

*4¹ means four seasons, Spring, Summer, Autumn and Winter. ² B is the abbreviation of blue bands of Landsat data, G for green bands, R for red bands, NIR for near infrared bands, SWIR1 for shortwave infrared 1 bands, SWIR2 for shortwave infrared 2 bands. ³ ASM is the abbreviation of angular second moment bands of textural features, CON for contrast bands, COR for correlation bands, VAR for variance bands, IDM for inverse difference moment bands, and SENT for sum entropy bands.

.

In terms of RF, there are two important parameters: the number of decision trees and prediction variables. The selection of these two parameters has a significant impact on the final results. The default number of decision trees is no less than the number of features of classification images, and the number of prediction variables is the basis of the feature number. These variables were set to 100 and 9, respectively.

3.3. Detection of Potential Area Changes

This step was designed to find potential area changes by calculating bi-temporal (basic year and comparison year) variables for each pixel rather than on an annual basis to avoid error accumulation. NDVI has proven to be less sensitive to varying sun–sensor geometry than other indices and is thus a robust metric for time-series analyses [48,49], and the change vector (CV) value has been frequently used to detect area changes, achieving satisfactory results [50,51]. Therefore, we calculated three variables: the CV value and the differences in NDVI maximum and NDVI minimum (dNDVImax and dNDVImin), as shown in Equations (1)–(3). The final potential area change was the union of the possible area change determined by these 3 indicators.

In this paper, $\mathrm{NDVImax}\_\mathrm{B}$ and NDVImin_B are the maximum and minimum NDVI composite images of the basic year; they are obtained by calculating the maximum and minimum NDVI values of each pixel. NDVImax_C and NDVImin_C are the maximum and minimum NDVI composite images of the comparison year (in Supplementary Materials code 2). For the CV, ${\mathrm{B}}_{Bi,max}$ and ${\mathrm{B}}_{Bi,min}$ are the reflectance of band i (blue, green, red, near infrared, shortwave infrared 1, and shortwave infrared 2) based on the maximum and minimum NDVI composite images of the basic year, respectively. ${\mathrm{B}}_{Ci,max}$ and ${\mathrm{B}}_{Ci,min}$ are the reflectance of band i based on the maximum and minimum NDVI composite images of the comparison year, respectively. Specifically, taking ${\mathrm{B}}_{Bi,max}$ as an example, we first determined the time at which the pixel obtained the maximum value of NDVI in the basic year. Then, we searched Landsat images according to the time. The band i value of this Landsat image at this pixel is ${\mathrm{B}}_{Bi,max}$.

$$\mathrm{dNDVImax}=\mathrm{NDVImax}\_\mathrm{B}-\mathrm{NDVImax}\_\mathrm{C}$$

(1)

$$\mathrm{dNDVImin}=\mathrm{NDVImin}\_\mathrm{B}-\mathrm{NDVImin}\_\mathrm{C}$$

(2)

$$\mathrm{CV}={{\displaystyle \sum }}_i{\left({\mathrm{B}}_{Bi,max}-{\mathrm{B}}_{Ci,max}\right)}^2+{\left({\mathrm{B}}_{Bi,min}-{\mathrm{B}}_{Ci,min}\right)}^2$$

(3)

The larger the values of dNDVImax, dNDVImin, and CV, the greater the possibility of change. Therefore, it is important to find the thresholds of dNDVImax, dNDVImin, and CV. We tested the mean plus 0.5 standard deviation, the mean plus 0.25 standard deviation, and the mean value plus 0.5 of the mean value as threshold values for the three indicators. The details of the results are shown in Figure 4. As we can see, new area changes or areas with slight differences may be excluded with thresholds set at the mean plus 0.5 standard deviation and the mean value plus 0.25 standard deviation. When the threshold value was 0.5 of the mean value, more unchanged regions were detected. The mean value was determined to be the most suitable threshold in this study.

3.4. Land Cover Mapping of Comparison Years

The outputs of the previous step were annual binary maps of the potential area change mask and unchanged area mask. The images of probable area changes can be obtained by using the mask of the potential area change to extract the image of the comparison year (in Supplementary Materials code 3). Then, we utilized RF to classify the images of possible area changes. The classification results were compared with the land cover types of the baseline year for each pixel. We updated the map of the baseline year with the results of the actual area changes and obtained the map of the comparison year.

3.5. Accuracy Assessment

The accuracy evaluation can be divided into two parts: the accuracy assessment of land cover classification results and multi-temporal change monitoring. For the first one, a confusion matrix was used to calculate users, producers, overall accuracy, and the kappa coefficient of each category. In order to test the accuracy of the change monitoring results, change samples and unchanged samples between 2010 and 2019, 2013 and 2019, and 2016 and 2019 were randomly selected at an interval gap of 3 years. High-quality samples were obtained based on high-definition images from Google Earth and Landsat data in GEE through visual interpretation.

3.6. Ecological Risk Assessment

ER reflects the impact of human activities on an ecosystem and can be reflected by the landscape pattern [29]. The ERI was established from the Landscape Disturbance Index (E_i) and Landscape Fragility Index (F_i) based on the land cover maps of 2010–2019. The calculation formula is shown in

Table 2. Construction method of Landscape Indices.

Name	Computing Method	Ecological Meaning
Landscape Fragmentation Index (Ci)	${\mathrm{C}}_i=n_i/A_i$	The higher the value, the higher the degree of fragmentation.
Landscape Isolation Index (Ni)	${\mathrm{N}}_i=\frac{A}{2A_i}\sqrt{\frac{n_i}{A}}$	The greater the value, the more complex the landscape distribution and the lower the ecological stability of the landscape.
Landscape Dominance Index (DOi)	${\mathrm{DO}}_i=\frac{Q_i+M_i}{4}+\frac{L_i}{2}$	The higher the value, the greater the impact of patches on the formation and evolution of landscape pattern.
Landscape Disturbance Index (Ei)	${\mathrm{E}}_i=a{C}_i+b{N}_i+cD{O}_i$	a, b and c are the weights of the corresponding landscape index, and the sum is 1. According to the existing research results and combined with the actual situation of the research area, the values are respectively 0.5, 0.3, and 0.2.
ERI	$ERI={\displaystyle {\displaystyle \sum }_i^N}\frac{S_{ki}}{S_k}\sqrt{E_i\ast F_i}$	ERI describes the size of ecological loss in the sample plot and converts the spatial pattern into ecological risk variable through sampling method.

n_i is the number of patches in the ith landscape; A_i is the area of the ith landscape; A is the landscape total area; Q_i is the sample number of the ith patches/ total number of samples; M_i is the number of ith patches/total number of patches; L_i is the area of the ith patches/ area of the sampels; S_ki is the area of landscape component I in unit k; S_k is the total area of unit k; and F_i is the Landscape Fragility Index.

.

E_i represents the degree of loss caused by the external disturbance of an ecosystem. This index is formed by the Landscape Fragmentation Index (C_i), Landscape Isolation Index (N_i), and Landscape Dominance Index (DO_i) [52]. C_i represents the degree of fragmentation of the landscape. N_i represents the degree of separation between different patch individuals in the landscape type. DO_i represents the important status of patches in the landscape. Landscape fragility is the vulnerability of ecosystem structure represented by various land cover classes. The fragility level can symbolize the ability of the landscape to maintain stability. The higher the score of F_i, the greater the ER. According to the characteristics of the study area and previous research [18,31], each type of land cover was given a corresponding weight. The weights of unused land, water and perennial snow area, cropland, forest land, grassland, and construction land were 0.2857, 0.2381, 0.1905, 0.0952, 0.1429, and 0.0476, respectively.

In landscape ecology, when the landscape sample (that is, the sampling unit) area is 2–5 times the mean area of landscape type patches, the landscape pattern information around the sampling point can be comprehensively represented [31]. Therefore, the scale of ER assessment was determined to be 3 km. The whole study area was divided into 3879 grids. The landscape index of each grid was calculated using Fragstats software.

Spatial autocorrelation was used to analyze high and low values of spatial aggregation [53,54]. Regions with high and low ER were identified with the Getis–Ord Gi* method. In addition, we determined the spatial distribution of ER by empirical Bayesian kriging.

4. Results

4.1. Multi-Temporal Classification Maps and Accuracy Assessment

4.1.1. Classification Results and Accuracy Assessment of the Basic Year

The overall accuracy of the 2019 land cover classification results was 93.67%, and the kappa coefficient was 0.92. The user and producer accuracies of water and perennial land were above 98%, as were those of cropland. The categories with lower recognition accuracy were grassland and unused land. The user accuracy of grassland was 83.69%, and the producer accuracy of unused land was 88.60% (as shown in

Table 3. Accuracy assessment for basic year classification.

		Ground Truth
		Unused Land	Water and Perennial Land	Cropland	Forest Land	Grassland	Construction Land	Users Accuracy (%)
Mapping Class	Unused land	443	0	0	0	0	31	93.46
	Water and perennial land	0	496	0	0	0	0	100.00
	Cropland	0	0	493	0	0	10	98.01
	Forest land	0	3	0	455	33	0	92.67
	Grassland	37	0	6	45	467	3	83.69
	Construction land	20	1	1	0	0	456	95.40
	Producer Accuracy (%)	88.60	99.20	98.60	91.00	93.40	91.20
	Overall Accuracy (%)	93.67		Kappa	0.92

). Commission errors mainly appeared between grassland and forest land and between unused land and construction land. In addition, some unused land pixels were misclassified as grassland. The detailed land cover maps are shown in Figure 5. From the perspective of distribution, most of the forest land and grassland are distributed in the south of Wuwei at high altitudes, the cropland and construction land mostly exist in the middle, and unused land is distributed in the north of Wuwei region. This is consistent with the basic situation of Wuwei region introduced in Section 2.1.

4.1.2. Classification Results and Accuracy Assessment of Comparison Years

Based on the classification results of the basic year, the changed regions from 2010 to 2018 were updated in the maps of comparison years. The accuracies of land cover maps in the comparison years were greater than 90%, and the kappa coefficients were above 0.88. In particular, in 2010 and 2018, the mapping accuracy was greater than 93%, and the kappa coefficients were above 0.93 (as shown in

Table 4. Accuracy assessment of multi-temporal classification maps.

Years	Overall Accuracy (%)	Kappa Coefficient
2010	93.83%	0.93
2011	90.27%	0.88
2012	91.63%	0.90
2013	91.10%	0.89
2014	92.40%	0.91
2015	92.30%	0.91
2016	92.33%	0.91
2017	91.40%	0.90
2018	93.83%	0.93
2019	93.67%	0.92

). The land cover mapping results of multiple years met the requirements for analyzing the land cover change pattern and ER. Classification errors mainly occurred between unused land and construction land and between forest land and grassland (the details are similar to those provided in Section 4.1.1) due to their similar regular spectral characteristics.

4.2. Area Change Detection and Accuracy Assessment

In order to investigate the effectiveness of the change monitoring method, a spatial random sampling strategy was used to select the changed and invariant samples. Intervals of 3 years, 6 years, and 9 years were used to evaluate the accuracy of change detection for 2016–2019, 2013–2019, and 2010–2019, respectively. The results for the three change monitoring periods were higher than 88.90%, which shows the effectiveness of the change monitoring method described here. The lowest accuracy of the change monitoring results was for between 2010 and 2019. The conversions between land cover types were complex, and the reflectance of the same category in 2010 and 2019 differed because the gap between 2010 and 2019 was bigger than that in the other two periods. Between 2010 and 2019, 17,287 pixels were randomly selected, including 8566 changed pixels. With the application of the change monitoring method in this paper, 7341 changed pixels were detected. The producer accuracy of the area change was 85.70%, and 8721 unchanged pixels were selected, of which 8072 pixels were detected. The producer accuracy of the unchanged region was 92.04%. The overall accuracies of change monitoring in 2013–2019 and 2016–2019 were higher than 93%, and the producer accuracies of the changed area and invariant area were higher than 91% and 95%, respectively. See

Table 5. Accuracy assessment of change area and non-change area.

		Ground Truth
		2010–2019			2013–2019			2016–2019
Mapped		Change	Non-Change	User Accuracy	Change	Non-Change	User Accuracy	Change	Non-Change	User Accuracy
	Change	7341	694	91.36%	9111	455	95.24%	5740	479	92.30%
	Non-Change	1225	8027	86.76%	866	9023	91.24%	564	10434	94.87%
	Producer Accuracy	85.70%	92.04%		91.32%	95.20%		91.05%	95.61%
	Overall Accuracy			88.90%			93.21%			93.94%

for details.

4.3. Land Cover Change Patterns

The areas of all land types were obtained from the land cover classification maps for 2010–2019 (shown in Figure 6). The area of each category varied greatly, and the most representative ones were unused land, cropland, and construction land. Overall, the areas of unused and construction land were reduced, while the area of cropland increased. Judging from the time series, the change in unused land area can be divided into two stages: 2010–2015 and 2015–2019. The area of unused land first decreased and then increased between 2010 and 2015, which also occurred between 2015 and 2019. The cropland and construction land areas showed opposite tendencies. In these 10 years, the area of water and perennial snow land increased significantly, from 122.67 km² in 2010 to 267.67 km² in 2019

To reveal the transition of land cover types in the study area from 2010 to 2019, the land use transition matrix was calculated. According to the results of the land use transition matrix (as shown in

Table 6. Land use transition matrix between 2010 and 2019 (Unit) km².

		2019
		Unused Land	Water and Perennial Land	Cropland	Forest Land	Grassland	Construction Land	Total
2010	Unused land	18,992.60	45.26	191.70	26.20	182.96	170.91	19,609.45
	Water and perennial land	2.69	117.90	0.03	0.19	1.73	0.14	122.68
	Cropland	78.80	11.62	3238.99	25.82	159.18	79.21	3593.62
	Forest land	22.43	46.13	2.97	2409.55	297.43	2.55	2781.06
	Grassland	65.99	44.20	17.41	229.26	2360.53	13.46	2730.85
	Construction land	102.47	2.56	395.74	25.27	58.13	3288.92	3873.10
	Total	19,264.99	267.67	3846.84	2716.10	3059.95	3555.20	32,710.75

), the main transitions occurred between construction land and unused land and between cropland and unused land. Unused land was mainly derived from construction land and cropland, with areas of 102.47 km² and 78.80 km², respectively. The degeneration of unused land was also significant, with areas of 191.70 km² and 170.91 km² converted into construction land and cropland, respectively. The most pronounced area change was in the water and perennial snow area. Although this area was not large, its acreage more than doubled from 2010 to 2019. The area of water and perennial snow was mainly converted from forest land, unused land, and grassland, with areas of 46.13 km², 45.26 km², and 44.20 km², accounting for 17.23%, 16.91%, and 16.51% of the water area and perennial snow area in 2019.

4.4. Spatial and Temporal Differentiation of Landscape Ecological Risk

Based on land cover classification products, we calculated the ER value of each pixel using landscape indexes (mentioned in Section 3.6). In this section, we analyze the spatiotemporal differentiation of risk from 2010 to 2019. The value of ERI was divided into five grades: lowest ER (ERI ≤ 0.209), lower ER (0.209 < ERI ≤ 0.231), moderate ER (0.231 < ERI ≤ 0.266), higher ER (0.266 < ERI ≤ 0.324), and highest ER (ERI > 0.324). The spatiotemporal differentiation of ER was significant in the study area.

From a spatial perspective, the distribution of ER in Wuwei showed a northwest to southeast spatial trend, which is similar to the distribution of topography and land cover types. From southwest to northeast, the distribution of ER can be divided into four zones: Zone1 (the lowest southern ER zone), Zone2 (the higher south-central ER zone), Zone3 (the lowest north-central ER zone), and Zone4 (the moderate northern ER zone), as shown in Figure 7a. The location of Zone1 is basically the same as that of the Qilian Mountain National Nature Reserve. During 2010–2019, the lowest ER area expanded, and the highest ER area decreased. Notably, the highest ER area was almost gone by 2019 (Figure 7). It was apparent that the area of the lowest ER had broadened, occupying almost the entirety of Zone1 (Figure 8). All signs in Zone1 indicated that the ER of the whole Qilian Mountain National Nature Reserve was gradually decreasing, and the ecological situation was improving. The area of unused and construction land accounted for the largest proportion of Zone2. With the implementation of EM, the areas of construction and unused land in this region decreased, while grassland and forest land expanded. Landscape fragility and fragmentation were relatively reduced. Therefore, the high ER area showed a disappearing trend. In Zone3, the connectivity of lowest ER area was enhanced, illustrating that low ER began to predominate. In addition, in Zone4, the area of the highest ER expanded, and that of the lowest ER shrunk, indicating that the ER value increased.

From the perspective of time, the average value of ER showed a downward trend from 0.2453 to 0.2398, with an overall decrease that reached 6.3%. This indicates that the ER of the study area was gradually reduced.

5. Discussion

5.1. Advantages and Limitations of Multi-Temporal Land Cover Mapping in GEE

With the help of open-source Landsat images and the powerful cloud computing capabilities of GEE, we detected changes in the EM area and obtained maps of land cover for 2010–2019. The open-source Landsat data laid a foundation for the study of multi-temporal change monitoring. A total of 2488 images were used in this experiment. Such a large amount of data is difficult to manually process and analyze on a computer with remote sensing software. The powerful cloud computing capability of GEE makes it possible to perform image preprocessing and image computing in large quantities. Additionally, there are many classification algorithm interfaces in GEE, which can classify a large number of remote sensing images with relative ease. Many advantages provide substantial convenience for the development of the remote sensing industry, but GEE has some shortcomings. First, for supervised classification, samples play an important role in model training and image classification results; however, when the number of samples exceeds 10,000, the memory limit is easily exceeded, which affects the accuracy of the classification results. Second, when using the cloud removal function in GEE, non-cloud pixels, such as snow, are more likely to be removed, resulting in unnecessary error. Finally, the pixel-based classification inevitably has noise points, but in GEE, the influence segmentation algorithm is not mature, and users cannot segment a wide range of images, which limits the use of the object-based classification method [55].

5.2. Analysis of the Method of Land Cover Mapping

The change monitoring method proposed in this paper was comprehensively improved from the aspects of characteristics and samples. To fully mine the hidden information of remote sensing images, the spectral, textural, terrain, and temporal characteristics were comprehensively considered. We randomly selected a determined number of sample points for each category from the training and validation datasets to ensure that the samples were representative and to avoid an imbalance in the classification results. In addition, potentially changed pixels rather than all pixels were classified so as to reduce the workload. Taking 2010 as an example, 61.63% of all pixels were identified as possible area changes, which needed further processing. Therefore, the workload of classification was greatly reduced, which highlights the efficiency of the algorithm in this study. False changes caused by classification errors could be significantly avoided, and the accuracy of unchanged categories was ensured. We opted for a more inclusive threshold in order to minimize potential omissions of any classes of land cover changes. Our results provide compelling evidence that this method can be used in other conditions of land cover change detection. When extracting the area change of the thematic map, other variables can be selected to find the possible area change, such as NDWI (in water extraction) or NDBI (in construction extraction).

5.3. Influence of EM

5.3.1. Influence of EM on Land Cover Change

From the analysis of the land cover change trend, the EM between 2010 and 2019 can be divided into two stages: the first stage (2010–2015) and the second stage (2015–2019). The Wuwei Municipal Government implemented the 12th and 13th Five-Year Plans in these two periods, respectively. The ecological resettlement project was an important part of each plan. At the beginning of each stage, a large amount of construction land and cropland was created to arrange for immigrations. In the development process of other cities and villages, construction land was mainly converted from cultivated land and green land [50]. The difference is that the construction land used for ecological resettlement in Wuwei was mainly derived from unused land. Another part of unused land was converted into cropland. Toward the end of each stage, a series of ecological construction measures was carried out. The cultivated and residential lands in the relocated area were abandoned or used for other purposes; however, the cultivated land and residential area were mostly transformed into unused land, and small areas of these categories were converted into grassland and forest land. This shows that the ecological restoration of the emigrated area has not been completed since EM, and it should be accelerated as a follow-up measure. The Qilian Mountains are an important functional water conservation area, which plays an important role in ecological regulation. The area of water and perennial snow increased. We believe that there are two main reasons. On the one hand, many reservoirs have been established to meet the needs of production and living because the study area is short of water resources. On the other hand, the rivers of the Qilian Mountains have been well protected during EM. The change in land types in the ecological migration area was drastic. In order to ensure the living and development needs of the residents, sufficient construction land and cropland were needed. At the same time, however, the green land as well as water and perennial snow area needed to be expanded to improve the ecological environment of the nature reserve. These requirements imply the significance of monitoring and managing land cover changes in the EM area.

5.3.2. Influence of EM on ER

The immigration and emigration caused changes in land cover, which, in turn, affected the dominance and fragmentation of the landscape. Additionally, different land types had different degrees of fragmentation. All of these aspects ultimately affected the ecological risks in the region.

In order to reveal the impact of ecological migration on the ecological environment of the study area, we selected the Qilian Mountain area and Huanghuatan resettlement site for analysis (Figure 9). The government has carried out a series of policy measures to restrict human activities in the Qilian Mountain area, contributing to the protection of the ecological environment. Due to a reduction in human activities, the areas of construction land and cropland have decreased. The fragmentation of the landscape has decreased. The predominance of grassland and forest land has been enhanced. The main landscape, circled in Figure 9d,e, has changed from unused land to water and perennial snow. This change greatly reduced the vulnerability of the area, benefited the restoration of regional habitats, and effectively reduced the ecological risk. It is undeniable that the EM project played an important role in water conservation in the Qilian Mountain area.

The Huanghuatan area is a vital resettlement site. To protect the livelihoods of the immigrants, large areas of residential land and cropland were built. In addition, a water storage and supply project was implemented. From 2010 to 2015, as the construction of the resettlement area advanced, large plots of land were divided into scattered areas of land; however, the infrastructure and supporting facilities were relatively complete. Therefore, the fragmentation of the landscape was reduced, while the stability increased. The ER of the resettlement area first increased and then decreased from 2010 to 2019, which is similar to the results of Liu et al. [9].

5.3.3. Influence of EM on Other Aspects

The results of this research indicate that EM has had a positive impact on protecting the ecological security of the study area; however, it might threaten the cultural diversity and ecosystem carrying capacity of the resettlement area. Many ethnic minorities reside in Wuwei as it is inhabited by 38 ethnic groups, e.g., Han, Tibetan, Hui, and Mongolian populations. After the implementation of the ecological migration policy, the production, lifestyle, ideology, and customs of ethnic minorities who originally lived in the impoverished area would have been changed, but the influence of these changes is difficult to judge. On the other hand, in order to guarantee the life satisfaction and well-being of migrants, a large number of infrastructure and supporting facilities need to be built in the resettlement area to ensure education, medical care, sanitation, transportation, employment, etc. This clearly reveals the general restrictive effect of the regional resource and environmental capacity on EM. Therefore, when formulating and implementing subsequent resettlement policies, the government must fully consider the resource and environmental capacity of the local area.

5.4. Limitations and Future Works

We constructed annual land cover maps from 2010 to 2019 using GEE and RF, and then we analyzed the impact of EM on changes in land cover and ER. There are still some uncertainties and limitations in this study. EM is a long-term process, and historical samples are difficult to collect through field surveys. Therefore, we selected samples through high spatial resolution images in Google Earth. It is inevitable that the results would have been affected by uncertain labels. In future research, we aim to investigate the study area and find more reliable official data to solve this problem. We selected ER as an indicator to quantify the ecological environment; however, the ERI describes the static pattern of ER rather than the dynamic process of risk adaptation and interaction [29]. It seems that using only ER is not sufficiently comprehensive. In further work, we plan to also consider the resilience and connectedness of the ecosystem. In addition, we will consider other factors, such as the ecosystem carrying capacity, ecological services, indicators of natural disasters, and risk sources of human disturbance in assessing the ecological environment. Establishing an evaluation system to quantify the impact of human activities on regional ecological risks is of great importance.

6. Conclusions

In this study, we developed a method for assessing the impact of ecological migration (EM) on changes in land cover and ecological risk (ER). Compared with other research focusing on human activities, our study has made significant improvements in ecological environment monitoring.

We constructed multi-temporal land cover maps of Wuwei for 2010–2019 with the help of the random forest algorithm. The characteristics and samples of the supervised classification method were comprehensively considered. Our approach entailed first monitoring and then classifying to reduce the extent of calculation. The accuracies of the land cover maps were above 90%, and the kappa coefficients were higher than 0.88 from 2010 to 2019. In order to further explore the impact of EM, we analyzed the transition of land cover types in the EM area. In addition, an ER assessment method was introduced to analyze the ER situation. The results show that: (1) the primary land use transition occurred between unused land and construction land and between unused land and cropland; and (2) the ecological condition of the Qilian Mountain area was protected by the EM project. During the EM project, the ER value of the ecologically protected Qilian mountain area declined. In the analysis of typical regions, we explored reasons to explain the ecological risk changes. This study provides a monitoring and evaluation method for the assessment of the effectiveness of EM on ecological environment. Our conclusions provide useful information and scientific guidance for efficient planning and sustainable development in the EM region. This idea can be extended to evaluate the impact of other similar human activities.