*2.3. Methods*

#### 2.3.1. MOLUSCE Plugin

Asia Air Survey Co., Ltd. (AAS) released MOLUSCE (Modules for Land Use Change Evaluation) at FOSS4G 2013, which was a conference for people working with open-source

tools. As a user-friendly plug-in for QGIS 2.0 and above, MOLUSE is designed to analyze, model, and simulate land use/cover changes. MOLUSCE is well suited to analyze land use and forest cover changes between different time periods, model land use/cover transition potential or areas at risk of deforestation, and simulate future land use and forest cover changes [63].

#### 2.3.2. Correlation Analysis

Correlation analysis refers to the analysis of two or more variable elements with correlation, to measure the closeness of the correlation between the variable factors. The measurement of the closeness of the relationship between geographical elements is mainly realized through the calculation and interpretation of the correlation coefficient. Pearson's correlation and Cramer's coefficient are the main correlation analysis methods in the MOLUSCE plugin of QGIS. Among them, Pearson correlation analysis is a measurement method of vector similarity [64]. The output range is from −1 to + 1, where 0 represents no correlation, negative value represents negative correlation, and positive value represents positive correlation.

The correlation degree is usually judged by the following value ranges:


#### 2.3.3. Change Analysis and Transition Potential Modeling

We used the MOLUSCE plugin inside QGIS to compute the land cover change between the research intervals. For transition potential modeling, we used the logistic regression approach. The elevation, relief, slope, monthly average temperature, annual average precipitation, river density, GDP, population count, road density, and city kernel density were used as the explanatory factors.

#### 2.3.4. Prediction and Model Validation

The MOLUSCE plugin can not only efficiently compute land cover change analyses, but is also well-suited for simulating future scenarios of land cover. We used the CA Simulation tool [65–67] of the MOLUSCE plugin inside QGIS to simulate the future land cover after we finished the transition potential modeling operation using the logistic regression approach. Next, we entered the reference map and simulated map for comparison and verification on the Validation page of the MOLUSCE plugin, and obtained the relevant Kappa coefficient values as a reference to check the accuracy of the simulation results.

#### 2.3.5. Annual Rate of Change Analysis

The annual rate of change (ARC) could be used to represent the magnitude of change between corresponding years. In order to obtain the annual rate of change for each land cover type, the area difference between the final year and initial year was divided by the area of initial year and time (year) period. We used Equation (1) to assess the annual rate of change in land cover categories [25,68]:

$$\text{ARC} = \frac{\text{Area}\_{\text{Final}} - \text{Area}\_{\text{Initial}}}{\text{Area}\_{\text{Initial}} \ast \text{t}} \times 100\% \tag{1}$$

where ARC is the annual rate of change in land cover categories. AreaFinal and AreaInitial are the areas of final and initial year, and t is the interval of years between the final year and initial year.

2.3.6. Landscape Pattern Index Analysis

The landscape pattern is the arrangemen<sup>t</sup> of landscape blocks of different sizes and shapes formed naturally or artificially in landscape space [69]. The landscape pattern indexes, as the sub-classification of landscape indexes, reflect the structural characteristics of land use/land cover types [70]. There are many types of landscape pattern indexes, and because of the application of new theories in landscape ecology, they are constantly being pushed forward [71,72]. Researchers often extend this part of the functions of the Geographic Information System (GIS) to form a unique landscape index software package based on GIS, such as Fragstats software package.

In order to study the spatial structure characteristics of different land cover types in the GYRR, this study first introduced the landscape diversity index to characterize them. Shannon's Diversity Index (SHDI) is a measurement index which is widely used in ecology based on information theory, and it is equal to the negative sum of the area ratio of each patch type multiplied by the natural logarithm of its value at the landscape level:

$$\text{SHDI} = -\sum\_{i=1}^{s} P\_i \ln P\_i \tag{2}$$

where s is the amount of patches, and Pi the area ratio of each patch type. When SHDI = 0, it indicates that the whole landscape is composed of only one patch, and an increase in SHDI indicates that the patch types increase or distribute equally in the landscape space. In a landscape system, the richer the land use/land cover is, the higher the degree of fragmentation is, and more uncertain information content leads to a higher calculated SDHI value. The diversity depends on two factors: the number of types and the evenness of area combination; therefore, the diversity index is the comprehensive embodiment of type richness and combination complexity [73].

Shannon's Evenness Index (SHEI) equals the SHDI divided by the maximum possible diversity under a given landscape abundance (all patch types are equally distributed). The smaller the SHEI value is, the more likely it is that some patch types may dominate the landscape, and a value that is close to 1 indicates that there is no obvious dominant type in the landscape while patch types are evenly distributed. Therefore, when SHEI = 0, it indicates that the landscape is composed of only one type of patch without diversity, and SHEI = 1 indicates that the patches are evenly distributed and have the greatest diversity.

$$\text{SHEI} = \frac{SHDI}{SHDI\_{\text{max}}} \tag{3}$$

where SHDI is Shannon's Diversity Index, and SHDImax is the maximum possible diversity under a given landscape abundance (all patch types are equally distributed) [74].

#### *2.4. Technology Roadmap*

In the process of this research, our work includes the following steps (Figure 5):

(1) We downloaded the land cover data for six years, including 1995, 2000, 2005, 2010, 2015, and 2020, and then, by comparing and analyzing the data of the first year (1995) and the last year (2020), we gained general insight into the land cover change in the GYRR in the past 25 years.

(2) The data of 10 geographical elements from different data sources, including elevation, relief, slope, annual average temperature, annual average precipitation, river network density, GDP, population, road density, and city density were collected.

(3) The MOLUSCE plugin was found in the plugin installation window and installed.

(4) We opened the MOLUSCE tool in the Raster menu drop-down list. Initial (2000) and final (2010) land cover data were used as input. Geographic impact factors, such as spatial variable, were used as input in the "inputs" tab of the MOLUSCE tool. Then, the subsequent operations were carried out step by step. The data generated in the previous

step were the basis for the next operation. In particular, we carried out the prediction of land cover in 2030 after verifying the effectiveness of the prediction model.

**Figure 5.** Technology roadmap. In this technical route, we not only paid attention to the land cover change in the whole GYRR, but also paid attention to the land cover change in mountainous areas which account for two-thirds of the whole GYRR, especially the change in residential areas in the mountainous areas.

(5) By comparing the land cover data in 2020 and 2030, we analyzed the land cover change over the next 10 years.

(6) On the other hand, we used the mountainous areas of the GYRR to cut out the land cover data of six periods from 1995 to 2020.

(7) The area changes in different land cover types in the mountainous areas of the GYRR were analyzed over time.

(8) Mountain hazard data points were downloaded and sorted to conduct the Euclidean distance analysis. The relationship analyses include distance from hazard points, number of settlement patches, distance from hazard points, and area of settlements.

## **3. Results**

#### *3.1. Correlation between Geographical Variables*

We calculated the Pearson correlation coefficient, as shown in Table 2. After comparison, it was found that the variables having strong correlation with each other include temperature and elevation, city density and temperature, city density and elevation, and relief and slope.


**Table 2.** Pearson correlation coefficient between different variables.

#### *3.2. Area Changes and Landscape Pattern Features*

The statistical analysis was done on various land cover areas between 1995 and 2020. The area change and the ARC of the same land cover were calculated (Table 3). It was noted that the land cover with the largest change was agricultural land, with a decrease of −16,437 km2. The increase in settlement area is the largest one, with an area of +27,364 km<sup>2</sup> and an ARC of 223.69%. The increase in settlement shows the enhancement of human activities in the past 25 years.

The area transfer analysis was also performed between different land cover types. According to the area transfer matrix (Table A1), between 1995 and 2020, large change situations include: 11,171 km<sup>2</sup> agricultural land was transformed into forest land, 53366 km<sup>2</sup> agricultural land was transformed into grassland, and 20,047 km<sup>2</sup> agricultural land was transformed into settlement land. In terms of forest land, 10,233 km<sup>2</sup> forest land was transformed into agricultural land, and 101,506 km<sup>2</sup> forest land was transformed into grassland. On the other hand, 52,088 km<sup>2</sup> grassland was transformed into agricultural land, 12,597 km<sup>2</sup> grassland was transformed into forest land, and 7645 km<sup>2</sup> grassland was transformed into settlement land. The Chord diagram (Figure 6) was used to express the land cover change. It was found that agriculture, grassland, and forest are the main land cover types, and account for most of the land studied.

We analyzed the landscape pattern indexes SHDI and SHEI in the GYRR, and the values of the two indexes increased gradually with time (Figure 7). The continuous increase in the SHDI value indicate that the land cover in the GYRR had become more and more abundant, and the higher the degree of fragmentation was, the greater the uncertain information content became. SHEI was getting bigger and bigger, approaching 1, which

indicated that there was no obvious dominant type in the GYRR, and landscape patches were more and more evenly distributed in the GYRR.


**Table 3.** Land cover change from 1995 to 2020.

**Figure 6.** Land cover change in GYRR from 1995 to 2020 using a Chord diagram expression. The right semicircle shows the proportions of different land covers in 1995, and the left semicircle shows the proportions of different land covers in 2020. The arrows in the circle indicate the land cover change.

**Figure 7.** Landscape pattern index analysis. (**a**) SHDI of GYRR; (**b**) SHEI of GYRR. The red line in the figure is a trend line added by the authors.

#### *3.3. Land Cover Prediction in 2020 and Validation*

We used the MOLUSCE plugin for the simulation of land cover in 2020. Using the projected 2020 data (Figure 8a) for comparison with the actual land cover data in 2020 (Figure 8b), the percentage of correctness was calculated as 96.42%, the Kappa (overall) was 0.94, the Kappa (histo) was 0.98, and the Kappa (loc) was 0.95. The results show that the 10-year interval prediction model has a good result on land cover simulation and prediction.

**Figure 8.** Actual and projected land cover in 2020. (**a**) Projected land cover in 2020; (**b**) actual land cover in 2020. We know that the more similar the two above maps are, the better the simulation results will be. However, there are still some subtle differences between the two maps. For example, the expansion trend of the simulated settlements was still conservative compared with that of the real settlements, and the real settlements expanded more rapidly, for example, in cities in Henan Province.

#### *3.4. Land Cover Prediction in 2030*

The above experimental results show that the 10-year interval land cover prediction model has good results. At the windows "Cellular Automata Simulation", the results show the option "Number of Simulation iterations". This means that, if only 1 is entered, it will be projected into the future only once. For example, if the land cover data are for 2000 and 2010, the land cover of 2020 will be projected when entering 1, and 2030 will be projected in the case of changing the "Number of Simulation iterations" to 2. In this study, the actual land cover in 2020 was used as the input of the model to simulate and predict the land cover in 2030 (Figure 9).

According to the statistical results of various land cover types, agricultural land, wetland, permanent snow and ice, shrubland, and sparse vegetation would be further reduced. The area of forest, grassland, settlement, bare area, and water would increase (Table 4).

**Figure 9.** Land cover prediction of the GYRR in 2030.

**Table 4.** Land cover change from 2020 to expected 2030.


In particular, it can be noted that a large amount of agricultural land would still turn into grassland and settlement. Some settlement land would be converted into agricultural land and grassland (Table A2). On the other hand, we analyzed the change in landscape pattern index (SHDI and SHEI) and found that the two indexes did not change much (Figure 10). This result shows that the land cover change in the GYRR may enter a stable development stage when it reaches a certain degree in the future.

**Figure 10.** Landscape pattern index (SHDI and SHEI) contrast analysis between 2020 and 2030.

#### *3.5. Land Cover Change in Mountainous Areas*

For the whole GYRR, the area percentage of plains and platforms is about 34.96%, and that of mountainous areas is 65.04% (Figure 2). The unique energy gradient causes the mountains to become an area of natural hazards development, such as debris flows,

landslides, collapses, avalanches, soil erosion, and mountain torrents. These mountain hazards may destroy urban and rural settlements, damage roads, bridges, and engineering facilities, bury farmlands and forests, and block rivers and reservoirs. They may cause huge casualties, property losses, and ecological damage, seriously threaten the lives and property of the people in the mountainous areas and the safety of engineering construction, and restrict the development of resources and economy in the mountainous areas [60]. We cropped out the land cover of the mountainous areas of the GYRR in different years (Figure 11).

**Figure 11.** Land cover in mountainous areas of the GYRR in different years.

Blind expansion of cities in mountainous areas can easily cause mountain hazards. It would cause huge losses to the life, property, and safety of urban residents. For example, on 14 August 2017, a devastating geo-hazard chain—debris slide, debris flow, and sedimentladen flood—occurred in Freetown, Sierra Leone, resulting in at least 500 deaths, more than 600 missing, and hundreds of houses destroyed. Although rainfall was a trigger factor for the Sierra Leone disaster, rapid and haphazard urbanization increased the hazard and vulnerability [75]. The development of mountain towns is generally affected by many factors such as social economy, topography, and geomorphology. Compared with plain towns, their infrastructure is relatively weak. In particular, poor urban planning and inadequate consideration of risks could lead to the construction of housing in dangerous areas. On the other hand, the removal of hillside vegetation increases erosion potential; low cost buildings using fragile building materials and methods could lack resilience; inadequate risk managemen<sup>t</sup> leads to weak emergency response.

We have also made area statistics for 10 land cover types (Figure 12). It can be seen that the settlement area in the mountainous areas had been increasing continuously in the past 25 years (Figure 12e), with an ARC value of +14.97%. Thanks to the policy of returning farmland to forests and grasslands implemented by the Chinese governmen<sup>t</sup>

in mountainous areas, the area of ecological land such as forest land and grassland had been continuously increased and, at the same time, the ecological environment had been improved as a gratifying result (Figure 12b,c).

**Figure 12.** Land cover in mountainous areas of the GYRR in different years. We tried to use colors similar to the colors of different land cover types.

We superimposed landslide, debris flow, and mountain torrent points on the base map of the GYRR (Figure 13). It was found that the mountain-hazard points are mainly distributed in the Central and Western regions of the GYRR. Next, the Euclidean distance is calculated using these mountain-hazard point data in order to represent the distance from the hazard point (Figure 13).

We made the statistics on the number of settlement patches, area of settlements in the mountainous areas of the GYRR, and the average Euclidean distance from the hazard points during the period from 1995 to 2020 with a time interval of five years (Figure 14). According to the statistical results, when the number of settlement patches in mountainous areas continued to increase, the distance between settlements in the GYRR and hazard

sites was also increased (Figure 14a). On the other hand, the mountainous residential areas in the GYRR also increased; however, the distance from the hazard point was also increasing (Figure 14b). The above two situations show that, although the intensity of human development in the mountains of the GYRR had been increasing, the awareness of avoiding hazards was also improved.

**Figure 13.** Mountain-hazard distribution and calculation of Euclidean distance around them.

**Figure 14.** Relationship between hazard points and settlements. (**a**) Distance from hazard points and number of settlement patches; (**b**) distance from hazard points and area of settlements.

## **4. Discussion**

*4.1. Decrease of Farmland and the Increase of Woodland under Returning Farmland to Forest and Grassland*

In 1998, Wuqi County of Yan'an, Shaanxi, China began to forbid grazing on the mountains [76]. After that, this small county, located in the northwest of the GYRR, began to take the lead in implementing the policy of returning farmland to forests [77]. Since 1999, Yan'an has reached a forest coverage rate of more than 50% and a vegetation coverage rate of more than 80% by returning farmland to forest over more than 20 years [78]. This is only a microcosm of China's project of returning farmland to forest and grassland. From 1999 to 2013, a total of 298,000 km<sup>2</sup> farmland in China was returned to forest. The project covers 2279 counties, with 32 million farmers and 124 million farmers directly benefiting. The central governmen<sup>t</sup> of China has invested 64.7 billion dollars in the project of returning farmland to forest [79]. From 2014 to 2018, the new round of returning farmland to forest and grassland involved 25 provinces (regions) including Hebei, Shanxi, Inner Mongolia, and others [80]. The Loess Plateau of the GYRR was one of the earliest regions to implement the project of returning farmland to forest and grassland which has made grea<sup>t</sup> contributions to the improvement of forest coverage in China [81]. Returning farmland to

forest and grassland provides a reasonable explanation for the reduction in farmland and the increase of forest land in the GYRR.

#### *4.2. Urbanization in Mountainous Areas of China*

Mountainous areas account for 24% of the earth's land area, and more than 12% of the world's population lives in mountainous areas [82,83]. China has a large number of mountains. The mountainous areas accounts for about 70% of the land area, and the population accounts for about 45% of the total population of the country. Due to the geographical and economic marginality of mountainous areas, the overall level of urbanization development in mountainous areas is far lower than the average level of China. The low level of urbanization and the slow urbanization process in mountainous areas, and a large number of agricultural population gathered in mountainous areas, will certainly bring grea<sup>t</sup> pressure to the ecological environment in mountainous areas. Among China's 1429 county-level administrative units in mountainous areas, 54% of the counties have an urbanization rate of less than 20%, and only 10% have an urbanization rate of more than 40% [84]. The development of natural resources, especially mineral resources, has played a significant role in promoting the urbanization of mountainous areas, and a number of resource-based cities have emerged. Tourism is a potential tool to promote the diversified development of mountain economy, increase the employment of mountainous residents, alleviate the poverty in mountainous areas, promote the participation of mountainous areas in economic globalization activities, and correct the development gap in mountainous areas. In recent years, tourism has become a new driving force to promote the urbanization of mountainous areas, such as Emeishan City, Wuyishan City, Tai'an City, and Jiuzhaigou County, and other counties and cities have developed rapidly through tourism. At the same time, due to the lack of managemen<sup>t</sup> and the lag in planning, some mountainous tourism cities have also experienced excessive urbanization [85]. Compared with plain towns, the urban planning in mountainous areas of China seriously lags behind the urban construction. At present, low-level spread of built-up areas and inefficient use of land are common in urban construction in mountainous areas. The level of urban functional layout is not clear. Especially with the increase in population and the shortage of construction land, the important functional layout of mountain towns basically ignores the avoidance of mountain hazards, and the ability of disaster prevention and mitigation is weak. For example, the area which was most seriously affected by the huge debris flow in Zhouqu, Gansu Province happened to be the most densely-populated and prosperous area [86,87]. The development of mountain towns in the GYRR also faces the above problems, accompanied by the increase in the area of residential areas and the number of residential patches.

#### *4.3. Harmonious and Sustainable Development of People and Land in Mountainous Areas*

The development of urbanization in mountainous areas should also be compatible with the resources and be coordinated with the land space and environmental capacity [88]. People should adhere to the fundamental support of ecological industry and form an intensive and ecological development model in order to improve the quality of urbanization. It is necessary to change the traditional direction and mode of urbanization development, gradually lead the urbanization construction to the road of new urbanization, implement the green development strategy, intensively utilize resources, improve resource efficiency, promote the intensive utilization of water, soil, and energy resources, and accelerate the construction of resource-saving cities and towns [89]. The development of urbanization in mountainous areas cannot ignore the restrictions of and close relationship with the mountainous environmental factors. We must grasp the basic characteristics and laws of the mountainous environment from different scales and regional differences, and deeply analyze the typical examples of the pattern, resources, and environmental characteristics, as well as the process of urbanization in mountainous areas. At present, the GYRR is carrying out the urbanization of the mountainous areas with an ARC value of +14.97%. There are many problems that need to be seriously considered to minimize the ecological

problems caused by the rapid urbanization of the mountainous areas and realize the harmonious and sustainable development of the relationship between people and land in the mountainous areas.

#### *4.4. Keeping Away from Hazards Benefitting from Disaster Prevention and Mitigation in Mountainous Areas of China*

China's mountains account for more than two-thirds of the total land area. With the rapid growth of population and inappropriate production activities in mountainous areas, mountain hazards occur frequently, and people have been disturbed and destroyed by mountain hazards in the process of utilization in mountainous areas [90]. Since the 1960s, China's relevant departments have begun to carry out the investigation and control of mountain hazards. For example, the color scientific and educational film "debris flow," released in 1965, brought the debris flow phenomenon onto the screen, which played a very strong role in publicizing and popularizing debris flow knowledge [91]. At present, many departments and colleges in China have trained many professional scientific and technological workers in theoretical research and disaster mitigation-prevention practice regarding mountain-hazards, and laid down a solid foundation in theoretical and technical reserves, becoming a very active scientific and technological force in China [92,93]. On the other hand, China has integrated the study of debris flows, landslides, floods, and other hazards, combined the construction of the large environment with the managemen<sup>t</sup> of small watersheds, carried out comprehensive research on the process of various hazards, implemented comprehensive disaster mitigation, scientifically assessed the current situation and trend of hazards, and put forward quantitative indicators [94]. With the increasing awareness of disaster prevention and mitigation in mountainous areas, although the proportion of residential areas and the number of residential patches in mountainous areas continue to increase, the safety of residential areas in mountainous areas has been continuously improved due to the conscious distance from mountain hazard points of urban construction [95].
