Monitoring Land Cover Change and Disturbance of the Mount Wutai World Cultural Landscape Heritage Protected Area, Based on Remote Sensing Time-Series Images from 1987 to 2018

Abstract

The contextual-based multi-source time-series remote sensing and proposed Comprehensive Heritage Area Threats Index (CHATI) index are used to analyze the spatiotemporal land use/land cover (LULC) and threats to the Mount Wutai World Heritage Area. The results show disturbances, such as forest coverage, vegetation conditions, mining area, and built-up area, in the research area changed dramatically. According to the CHATI, although different disturbances have positive or negative influences on environment, as an integrated system it kept stable from 1987 to 2018. Finally, this research uses linear regression and the F-test to mark the remarkable spatial-temporal variation. In consequence, the threats on Mount Wutai be addressed from the macro level and the micro level. Although there still have some drawbacks, the effectiveness of threat identification has been tested using field validation and the results are a reliable tool to raise the public awareness of WHA protection and governance.

1. Introduction

A world heritage site is a site with special historic, scientific, or esthetic qualities, and a universal value [1]. This has been widely accepted since the concept was first raised by the United Nations Educational, Scientific and Culture Organization (UNESCO) in 1972 [2]. There are four kinds of world heritage, as follows: Natural heritage, cultural heritage, mixed heritage, and cultural landscape heritage [3]. While natural and mixed heritage are closely related to the environment, cultural heritage seems more independent [4]. Cultural heritage includes tangible cultural heritage, such as movable, immovable, and underwater cultural heritage and intangible cultural heritage, which are treasures of human activities, for instance, oral traditions, performing arts, and rituals [5]. These kinds of cultural heritage interact deeply with the natural environment [6]. It is meaningful to place the cultural heritage within the surrounding landscape together, as the pattern and context of this system will influence the function and integrity of the whole ecosystem directly [7]. Thus, the analysis of the World Heritage Area (WHA) should encompass heritage sites and their surrounding landscape.

However, environmental change has posed a great disturbance, and even threats to WHA’s security, and has raised a number of problems concerning its sustainable protection in a specific area. For instance, excessive tourism may destroy the balance of the WHA [8]. Continued urbanization not only causes land resource scarcity, but also leads to an urban heat island [9], which directly changes the characteristic of the heritage [10]. Large-scale agriculture, mining activity, and erosion could also threaten the development of the WHA [11]. In other words, monitoring WHA and their surrounding environment at the regional level is a fundamental part of cultural heritage protection.

Moreover, an ecological disturbance can be described as “a case; a physical force, agent, or process, either abiotic or biotic, causing a perturbation (which included stress) in an ecological component or system; relative to specified reference state and system; defined by specific characteristics” [12,13]. As kinds of protected areas, WHA have shown common features and own unique characteristics [12]. It should be noted that both WHA and protected areas have suffered serious ecological degradation, which threatens the balance of WHA and the whole ecological system badly [8,14]. For instance, species’ population and habitat conditions have decreased sharply in the Great Barrier Reef Marine Park due to the poor management [8]; LULC change, caused by urbanization or mining, generates a significant influence on WHA [15]; the Arabian Oryx Sanctuary in Oman has been deleted from the world heritage list because oil and gas extraction have occupied almost ninety percentage of its original area [8,16]. The land pressure by urbanization and mining activity leads to complex problems in the ecosystem [17,18,19]. Meanwhile, social, economic, scientific, and political changes in the protected areas are closely connected to the processes of protection [15,20]. In some developing countries, tourism is the main source of income for the development of the local economy. For instance, Rwanda earns about 200 million US dollars per year from tourists visiting the Volcanos National Park [8,21].

No matter how pollution derived parameters or land use/land cover (LULC) change, the impact of constraint factors on the WHA is a prolonged process [12,22]. On the one hand, there is evidence that human activities, such as mining, cutting trees, and fossil fuel overuse significantly increase the generation of acid gases and other air pollutants [23,24,25]. It is obvious that, under such an environmental background, the WHA will suffer accumulating damage in the long run [6]. On the other hand, LULC change means degradation or increase for a specific category and it will produce a direct spatial pressure on WHA [26,27]. Population migration, technological development, economic growth, politics, and values are the driving forces of LULC change [11,25]. According to studies based on the pressure state response indicator framework [27], LULC change has a significant influence on the regional water cycle, environmental quality, biodiversity, productivity, and adaptability of the terrestrial ecosystem. Meanwhile, the density of the built-up area and vegetation cover have also vastly impacted land eco-security [28,29,30,31,32]. It is certain that the LULCC is a continued and dynamic process. Therefore, exploring the impact of LULC in WHA from a time series perspective is quite useful to comprehensively understand the process.

Remote sensing (RS) techniques have revolutionized the traditional recording and prospection approach [7]. By adopting such well-established techniques, the measurement of geographic units becomes more efficient and detailed [33]. Firstly, RS techniques could provide significant data support for environmental assessment, planning, modelling, and prediction [34]. For instance, RS image classification could use medium or high spatial resolution image data and various classifiers, including support vector machine, rotation forest and multi-classifier, to obtain the LULC information [35,36,37]. In order to achieve a comprehensive analysis of the dynamic changes of the WHA environment and its driving factors, time series RS products could illustrate the changing process and mechanism [38]. Furthermore, RS techniques could calculate and simulate plenty of key environmental indicators, including chlorophyll a [39], suspended solids concentration [40,41], surface radiance energy budget [42,43], surface moisture and temperature changes [44,45], vegetation cover and land degradation [46,47], vegetation biomass [48,49], net primary production (NNP) [50,51], vegetation structure [52,53], ecological parameters [51,54], surface ecological process [55], land use and agricultural vegetation extraction, and land cover vegetation ecological process [56,57,58,59,60,61].

In terms of the RS assessment systems of protected areas, there are several frameworks based on different methods. Firstly, with the system based on PCA (principle component analysis), for instance, RSEI (remote sensing based ecological index) determined greenness, dryness, wetness and heat, four main ecological indicators, then (PCA) was utilized to compress the four indicators into one in order to assess overall regional ecological status [62]. Secondly, there is a system based on weights, which come from experiments or experts, such as applications in groundwater potential zones [63], human settlement environment [18,64,65], mix research in environment and society [66,67,68], etc. Then, there is the system based on pattern recognition methods, such as hierarchy analysis (HA), analytic hierarchy process (AHP), fuzzy evaluation (FE), genetic algorithm (GA), and statistical learning (SL) [65,69,70,71,72].

In conclusion, RS assessment systems could be an efficient method to monitor protected areas. According to Reference [11], research objective has been divided into four categories, as follows: Habit mapping and change detection, the assessment of the habitat degradation, the evaluation of species’ diversity and distribution, and the tracking of pressures and threats. Up to now, there are a large number of RS assessment systems that have been established to monitor protected areas [58,73,74,75]. There are several assessment systems for different purposes, however, they follow a similar structure, which consists initially of detecting disturbances in the protected areas. These disturbances can be disturbances in identification, extraction indicators, and assessment. For instance, firstly, disturbances identified in protected areas generally include urbanization, roads, construction, mining, logging, agriculture, fire, alien species, hunting, grazing, and drought. Secondly, extraction indicators derived from remote sensing datasets have expressed disturbances to the WHA or the eco-system. Then, the assessment framework can be constructed to measure the pressures and threats that influence the protected areas. Although most WHAs have archaeological value, there is no doubt that the protection of WHAs is imperative and urgent [76].

This paper focuses on one specific WHA and seeks to identify pressures and disturbances within it, aiming to better understand the interaction between WHA and the surrounding environments. To conduct the research, a comprehensive assessment framework for the WHA is proposed, which analyzes the pressures and threats from every single aspect and overview perspective. Based on human activities and the natural environment of the selected WHA, four potential disturbances (forest coverage, vegetation condition, built-up area, and mining area) are mainly considered in the assessment. To assess the WHA at a regional level, this paper has adopted a synthetic pressure assessment method which integrates three key indicators (land stress, vegetation coverage, and water network denseness) and the biological richness derived from LULC. Meanwhile, the time series analysis was successfully employed for describing and mapping the spatiotemporal evolution of disturbances and WHA threats from 1987 to 2018. As the disturbances of WHA have differences in spatiotemporal levels, to comprehensively identify the characteristics of spatiotemporal variation and distribution, the linear regression analysis combined F-test method was adopted. There are serval studies that focus on the Mount Wutai environment, but few of them were based on a time-series remote sensing technique. This research monitors both LULCC and disturbances around Mount Wutai by adopting a remote sensing technique and proposes a unique index (CHATI) to comprehensively assess the threats in the research area. In addition, as we mentioned that a lot of RS assessment systems have been employed in environment assessment, however, few of them only use RS indicators, which reduce the data dependency in the accepted error range. Furthermore, in this research, the spatiotemporal analysis is based on linear regression and F-test, which provides more detailed information of spatial variation.

2. Study Area and Data 2.1. Study Area

Mount Wutai, located in the southeastern Shanxi Province, China, has a very unique geological structure, magnificent natural scenery, and abundant cultural significance. It mainly consists of two parts, Taihuai Proposed Core Area (N 39°01′50″, E 113°33′48″) and Foguang Proposed Core Area (N 38°52′56″, E 113°20′58″), and the core area has reached 18,415 hectares, which is composed of a series of mountains. The lowest elevation of Wutai Mountain is only 624 meters, while the Beitai Ding, known as “the Roof of North China”, has an altitude of 3,061m above sea level. The unique geological structures, stratigraphic profiles, fossil relics of paleontology, Cenozoic planation, and periglacial landforms from the early stage of the earth are perfectly preserved in Mount Wutai. On the mountain, the temperature is relatively low, such that snow covers the high peaks for years. So, it is hard to identity the landscape in winter by RS. The natural resources are well protected and 8.66% of the mountain is covered with forests with vertical pines, firs, poplars, willow trees, and lush grassland. The pleasant climate and beautiful scenery attract countless tourists and artists [77,78,79]. As one of the four sacred Buddhist mountains in China, it was affirmed as a UNESCO World Landscape Heritage Site in 2009. It is considered as the bodhimanda of Manjusri (Wenshu in Chinese), who is the Bodhisattva of wisdom. It is said that Buddhism was introduced to Mount Wutai in the Eastern Han Dynasty. Since then, it has become one of the most important centers for the development of Buddhism and has made significant contributions to the development of the Buddhist culture [77,78,79,80,81,82].

With the development of the economy, human activities play more and more important roles in local environment. As Mount Wutai has abundant tourism and mineral resources, the development of tourism and the mineral industry has become the focus of local economy. Here, take the development of tourism as the example. Since 1982, Mount Wutai has been affirmed as national scenic spot by the State Council and the local government adopted a series of measures to develop tourism, which included strengthening the construction of tourism infrastructure and building a relative tourism supporting service system that integrates food, housing, transportation, travel, shopping and entertainment [83]. The policy of prioritizing tourism had achieved remarkable economic benefits. Until 2008, a total of 2.8102 million domestic and foreign tourists were received throughout the whole year, and the total income of tourism reached 1.405 billion yuan [84]. Over the past three years, the number of tourist receptions continued to rise, with an average of more than 4 million per year. It can be said that the revenue from tourism and the mineral industry has already supported half of the local government’s financial benefits [85]. However, such tourism activities, the construction of tourism infrastructure, and developments in the mineral industry have already seriously affected the local ecological environment, such as the variation of the forest coverage. The influences of human activities on the study area are analyzed in detail by adopting remote sensing techniques.

UNESCO demarcates Mount Wutai into a core zone and a buffer zone. The Taihuai Proposed Core Zone refers to the area centered around the temple ensemble at Taihuai Town and the five terrace tops of Mount Wutai. The Proposed Buffer Zone is the protected areas of Mount Wutai. To identify the landscape and threats to Mount Wutai, the research area is demarcated by administration boundary and mountain ridges (as shown in Figure 1). The total area of the research area is about 4781.8km².

2.2. Data

This research uses multi-source and multi-temporal data. Firstly, MODIS NDVI (normalized difference vegetation index) products (MOD13Q1), generated by NASA, are used to identify the seasonal variation of vegetation coverage and growth situation. A total of 431 phase MODIS NDVI product images, from 2000 to 2018, with 250 meter spatial resolution and 16 day temporal resolution were used to derive the mean NDVI value over Mount Wutai (Figure 2). Such results provide guidance to the image selection of Landsat data.

Secondly, the main data in this research are Landsat time-series data. Those data captured from 1987 to 2018 by Landsat sensors, mainly including TM (Thematic Mapper) and OLI (Operational Land Imager) images with 30 meter spatial resolutions and less than 90% cloud cover (based on a visual estimate). According to the mean NDVI of Mount Wutai from 2000 to 2018 (Figure 2), resulting from MODIS NDVI products, regular patterns of seasonal dynamic are shown in the research area. Thus, the data during the vegetation growing season from June to September are considered as the best choice to observe the vegetation situation. However, for some years, such as 1995, 1998, and 2018, we could not obtain cloud-free Landsat data, we choose adjacent month’s data. In addition, to show the changes in significance, we choose data every 2 years. Overall, a total of 12 phase data are shown in Table 1.

Thirdly, slope data and aspect data generated by ASTER DEM with 30 meter spatial resolutions are used as ancillary data to improve the classification accuracy.

Finally, 1:100,000 topographic map data and forest resources inventory data were collected from 2016 to 2017 and were used in the field survey.

All TM/OLI data are processed by ENVI and ArcGIS software. During the image preprocessing period, radiometric calibration is used to eliminate the error of the sensor itself and to determine the exact radiation value at the sensor, and atmospheric correction is adopted to eliminate errors caused by atmospheric scattering, absorption, and reflection. Geometric correction is employed to eliminate the effects of terrain or deformation caused by camera orientation, mosaic, and subset.

3. Method

Long-term time series data have proven to be useful in detecting LULCC, and the subsequent environmental modelling across the regional and global scales [86,87,88]. Some previous studies have successfully used time-series data in the protected areas [11,89,90,91]. In the research, 12 phases of TM/OLI from 1987 to 2018 are the main time series dataset. Based on the indexes (NDVI, NDSI (normalized soil index), MNDWI (modified normalized difference water index)) and LULC data derived from Landsat time-series data in Mount Wutai, this paper identifies the potential risks and major disturbances of the WHA, including biological abundance, land stress degree, vegetation coverage, and water network density, and the threat is assessed comprehensively. Firstly, GIS and RS techniques are combined with multi-source data to extract those key factors from Landsat time-series imagery. Secondly, the typical characteristics of the heritage site are analyzed from a single factor perspective, then a comprehensive evaluation system and the long-term sequence framework are established. The comprehensive risk of the heritage site is further evaluated. Finally, the result of single factors, such as the vegetation, and a comprehensive evaluation are carried out across the perspective of spatial and temporal scales and the changes in the heritage site are analyzed throughout the micro and macro scales. The overall method is shown in Figure 3.

3.1. LULC and Indices Extraction

The LULC and indices (NDVI, NDSI, and MNDWI) are the fundamental baseline data of this research. They are partially the key elements of land cover disturbances identification and all contribute to the comprehensive assessment framework. As shown in Figure 3, the LULC is derived from multi-source data, which includes the following steps: (1) The selection of Landsat data based on seasonal variation of vegetation coverage and growth situation, with statistics from MODIS NDVI products in Mount Wutai; (2) making a radiometric correction, an atmosphere correction, a geometric correction, and a subset; (3) combining the ASTER DEM, slope data, the and aspect data with all bands of processed Landsat time-series data; 4) using 3 different single classifiers, maximum likelihood classification (MLC), neural net classification (NN) and support vector machine classification (SVM), to achieve 3 LULC classification results for each selected year. The classifications are able to differentiate 7 LULC categories, as follows: Forest, grass, farmland, building, mining, soil, and water. Based on validation, the SVM, MLC, and NN show good accuracy in green space, buildings, and bared soil in the study area. The 3 LULC classification results are fed into a classifier ensemble, where a simple majority or weight voting is chosen as the ensemble method. The ensemble classifier is the most suitable approach to capture multi-source and complementary information within each individual classifier, which is the same principle as with Random Forest classifier, with respect to single decision trees [35]. A classification accuracy assessment is made to validate the results rigorously.

The main indices used in this paper are NDVI, NDSI, and MNDWI. NDVI (Normalized difference vegetation index) is the ratio parameter of the near-infrared band (NIR) and the infrared band (R) reflectance of remote sensing images. The MNDWI (Modified Normalized Differences Water Index) is developed from the NDWI [46], which can effectively extract water information from remote sensing images. The NDSI (normalized soil index) is developed from the Soil Index (SI), where SI (Soil index) is combined with IBI (index-based building index), and uses SAVI (soil-adjusted vegetation index) to exclude the vegetation influence. Consequently, NDSI is the suitable indicator for the building area and bared soil. The comprehensive detail of the indices (equations and further explanations) are listed in Table 2.

3.2. Identification of Land Cover Disturbances

Disturbances can be divided into a variety of categories. In general, satellite-based remote sensing focuses on urbanization, road construction, mining, agriculture, fire, invasive species, hunting, grazing, and drought. For the Mount Wutai Heritage area, as mentioned in the introduction, the main disturbances concerned forest coverage, built-up area, mining area, and vegetation conditions. The first three disturbances are derived from LULC data, and the vegetation conditions are characterized by NDVI [53]. Regional statistical approaches are employed to quantify the forest coverage rate, the built-up area, and the mining area. The vegetation index is used to characterize the growth of vegetation in the region. According to Zhu (2013) [100], four vegetation indices, the normalized difference vegetation index (NDVI), simple ratio vegetation index (SR), the modified normalized vegetation index (MNDVI), and the reduced simpler ratio of vegetation index (RSR), are compared in mountainous forests using Landsat TM data. Their research suggests that terrain is strongly influenced by ratio vegetation index (RSR and MNDVI). Meanwhile, SR and NDVI can largely eliminate the influence of terrain. NDVI shows a better performance in this area. Therefore, NDVI is selected to characterize the vegetation growth in this study area.

3.3. Comprehensive Assessment System of Threats in the Heritage Area.

The comprehensive assessment system of threats in the heritage area is developed from the EI (Environment Index), which is constructed according to the “Technical Specifications for Evaluation of Ecological Environment Conditions” of the Ministry of Resources and Environment, China [101]. In this system, ecological environment conditions can be evaluated from six aspects, including the biological richness, vegetation coverage, water network denseness, land stress, pollution load, and environmental restriction. To systemically assess the environment conditions, each of the six aspects developed an index with a corresponding weight. This weight has been published as a guide book, as shown in Table 3.

In terms of the remote sensing assessment, the EECI (Ecological Environment Conditions Index) has been developed.

EECI=α×NHQ+β×NNDVI+γ×NWD+δ×NLS.

This system is focused on the measurable indicators; thus, the environmental restriction index will not be considered here. In addition, the pollution load index is characterized by chemical oxygen demand, ammonia nitrogen, sulfur dioxide, dust, nitrogen oxides, and solid waste in the EI. These indicators show the bearing capacity of the environment for a specific area. |Since these data cannot be effectively captured by remote sensing techniques at some specific areas, the bearing capacities always keep stable, despite some emergency situations. The weight of the pollution load index can be a close to 0. Moreover, in the protected areas, vegetation situation could be more sensitive than in other areas. Thus, the weight of vegetation coverage index can trend to 0.35. Finally, the EECI consists of the biological richness index, the vegetation coverage index, the water network denseness index, and the land stress index. Overall, the coefficients in this equation are similar to EI, α is 0.35, β is 0.35, γ is 0.15, and δ is 0.15.

The NHQ (normalized habitat quality index) is characterized by spatial data after bio-abundance, using LULC type reclassification. NNDVI (normalized NDVI) is characterized by normalized vegetation index (NDVI). NDVI is extremely sensitive to the surface vegetation and vegetation growth situation. Thus, it is considered an effective indicator for monitoring regional vegetation and ecological environment changes. The MNDWI has been tested on remote sensing images with different water types, and most of them achieved better results than NDWI, especially the extraction of water bodies within the urban area. Additionally, MNDWI can easily distinguish between shadows and water bodies, which solves the difficulty of eliminating shadows in the water extraction problem. NWD (normalized water network denseness) can be characterized by MNDWI. NLS (the degree of land stress) is the land stress to the environment at the regional level. The stress comes not only from bared soil, but also from the built-up area. NDSI combined IBI with SI can avoid the loss of information on construction area by applying SI, in such a case that NLS is characterized by NDSI.

In terms of the EECI, a high EECI value indicates the ecological environment conditions are more sustainable and less influenced by the disturbance. Instead, a low value means a high threat in this area. In order to analyze the result consistent with the threat value, the CHATI has been defined as follows:

CHATI=1−(α×NHQ+β×NNDVI+γ×NWD+δ×NLS).

The WHA threat assessment index is the indicator to determine the existing risks of the WHA at a specific time. The value of CHATI is calculated from a lot of spatial variables. The weights for each parameter and indicator come from expert experiences and domain knowledge. According to the weights of EECI and multiple tests, the final weights are shown in Table 4.

3.4. Variation Trends of NDVI and CHATI

As the ecological environmental comprehensive assessment indicators are spatially different, in order to systematically explore the spatial variability and space-time characteristics of the comprehensive ecological environment indicators in the regional level, the method of linear regression analysis [102] + F test [103] is wildly used in spatial analysis and can provide reliable results [104,105,106], with the year as the independent variable. To determine the eco-environmental quality of the CHATI and NDVI indices as dependent variables, a linear regression equation for each pixel is constructed in the study area as follows:

_θslope=n×∑i=1ni×_Ci−∑i=1ni∑i=1n_Cin×∑i=1nⁱ²−^(∑i=1ni)2,

where θ_slope is the slope of the regression equation for NDVI and CHATI; n is the number of years of the study period (11a for 1987–2018); i is the year number from 1 to 12; and C_i is the NDVI and CHATI sequence data for the study subjects. If θ_slope is positive, it indicates that the research object changes with time, and the larger the value, the more obvious the upward trend. Otherwise, the research object declines. The significance test of the trend can also be tested. The significance only represents the level of confidence in the trend change and has nothing to do with the seed of change. Its calculation formula is as follows:

F=U×n−2Q,

U=∑i=1n(y^i−y¯)²,

Q=∑i=1n(yi−y^i)².

In the formula, the U is sum of the squares of the regression; Q is the sum of the squares of the errors, in which y_i is the actual observation value of the year; i andŷ_i is the regression value, which is the average value of 12 years; whileӯ_i is the number of years. According to the test results, the trend is divided into the following seven levels: Significant decrease (Slope < 0, p < 0.05); decrease (Slope < 0, 0.05 < p < 0.1); slightly decrease (Slope < 0, 0.1 < p); slightly increase (Slope > 0, 0.1 < p); increase (Slope < 0, 0.05 < p <0.1); and significant increase (Slope < 0, p < 0.05)

4. Result and Analysis 4.1. The Distribution and Variation of Main Disturbances

Land use/land cover change-based remotely sensed time-series imagery has been widely used in land resource monitoring and management [107]. This method can select different classification categories depending on the research objectives [108]. In this research, we focused on the forest, grass land, mining area, farmland, building, water, and bared soil to understand the disturbances in WHA. The LULC maps from 1987 to 2018 are generated by a supervised classifier ensemble, as shown in Figure 4.

Those LULC maps obtained from remote sensing data are accompanied by a certain error probability [109]. There are many reasons that could lead to error probability, such as the selected classifier (based different algorithms), image quality (cloud coverage, sensors type and so on), and training data quality [107]. To assess the classification results, the confusion matrix including overall accuracy (OA) and Kappa coefficient [110] were adopted in this research (Table 5). From this table, the majority of the OA is over 85%, which is considered as the acceptable classification accuracy in remote sensing applications [111].

To better understand the changes between all categories, 11 from-to tables that include the area of the LULC classes for each selected year have been listed (Table 6).

From the data statistic (Figure 5) and LULC maps, the development of Mount Wutai over 30 years, from 1987 to 2018, is demonstrated comprehensively. Overall, forest coverage of Mount Wutai is in decline. On the contrary, building area shows a significant increase. The mining area exploded to a peak point and then slowly decreased. The bared soil shows an increasing trend, whereas the farmland decreased initially, with a slight increase thereafter. A small area of water is found in Mount Wutai, where the area keeps stable. In terms of the grassland, the area fluctuates drastically. Detailed descriptions of each category and disturbances are illustrated as follows.

Firstly, this paper focuses on the four main disturbances of this area which has been identified in the methodology, including forest coverage, vegetation conditions, built-up area, and mining area, respectively.

Forest coverage is derived from the LULC map. From the macro-perspective, the forest coverage of Mount Wutai has declined in fluctuation. In 1982, Mount Wutai was approved by the State Council to be included in the list of the first batch of national scenic spots [81]. It raised the public’s awareness on protecting green vegetation. At the same time, Wutai County was only 0.2 billion RMB Yuan in 1987, according to the GDP references (Gross Domestic Product) [112]. The secondary industry only occupied 14.5% [83], where the majority of people lived in the countryside [112]. Thus, the human activities had little influence on the environment of Mount Wutai. The forest coverage remained stable at that time. In 1992, Mount Wutai became a national forest park. The government of Wutai County received the special funds to plant more trees and protect the environmental of Mount Wutai [81], where forest coverage reached a peak point. However, since the 1990s, industrialization has increased across the whole Shanxi Province, China [113] (Figure 6). Human activities became the main factor of the regional and environmental balance [85], as shown in the from-to Table 6., while some of the forests changed into a built-up areas or mining areas. In 2000, the government of Mount Wutai prepared to apply to be a World Heritage Site [114]. A tree planting program began to show a positive impact of forest coverage. In only 2017, 3,314,863 hectares of trees were planted in the Mount Wutai area [113]. From the micro perspective, the forest in the LULC maps shows an increase around the core protected area of Mount Wutai. The main degradation area is concentrated in built-up areas and mining areas. It is consistent with the existing literatures [83,112,113,114].

For the vegetation conditions, there is a large number of remote sensing indices that have the ability to efficiently and effectively characterize vegetation conditions [46,93,99,100,115]. In this paper, NDVIs generated from Landsat images are adopted for each selected year to indicate the vegetation conditions. For the entire Mount Wutai, the mean NDVI kept stable, despite the seasonal effect in 1995, 1998, and 2018 (Figure 7). To identify the spatial variation of vegetation conditions, a significance of variation in NDVI, based on an F-test, is shown in Figure 8. This figure can directly show the trend in changes of each pixel. As the figure shows, the vegetation condition significance decreases in the mining area and built-up area, and the slight red line in the map is consistent with the road network that has been built over the last 30 years.

Meanwhile, the built-up area fluctuated dramatically, as shown in Figure 9. Over the period from 1987 to 2001, the number of households in the Shanxi province increased at a fast pace [84] (Figure 10). However, the urbanization process has not significantly influenced Wutai county, due to the slow economic development, and the area is stable between 50 Km² to 100 Km². In 1994, the central government of China made a decision to reform the urban housing system [116]. This event transferred the urban house system from public housing system into housing market system [117]. After that, crowed real estate companies rushed into the province-level house market and the county-level house market [118]. The years 2001–2004 witnessed a sharp rise in the built-up area and, in 2004, the built-up area became more than 170 Km². Another reason for which is because fundamental infrastructure investment sharply increased in the whole Shanxi Province [113] (Figure 11). For Mount Wutai, the total road length increased from 100 Km in 1987 to 1227 Km in 2017 [119]. However, after 2004, the population of the Mount Wutai area increased slightly [120], due to an increasing number of people migrating to the large cities [121]. In addition, the investment into fundamental infrastructure was also abandoned. The built-up area began returning to a stable stage.

Mount Wutai has abundant mineral resources, including 30 kinds of metal and non-metal mineral products, among them, iron ore and coal are the most common [112]. Consequently, the exploitation of mineral resources is also an important part of the development in Mount Wutai. From 2011 to 2014, more than half the tax was from the industry [120] (Table 7). Owing to the limitation of the economy and technology, the mining area remained steady from 1987 to 2001 and has not surpassed 50 Km².

In the year 2001, due to the high demand for coal and iron, all of Shanxi Province expanded the production capacity for industrial production [113] (Figure 12) and, consequently, the economy of Wutai County developed rapidly [119] (Figure 13), with a fundamental shift towards growing energy markets. The mining area in Mount Wutai expanded dramatically, which is the second reason for the downturn of the forest coverage, and reached a peak at almost 300 Km² in 2007.

During that time, a lot of illegal surface coal mining appeared [122,123]. In 2008, the energy market of Shanxi Province has suffered from a cold winter, due to the world economic crisis [124,125]. The development of mining areas tends to be stable. During the following three years, the mining area has gone through a small amount of decline. From 2010 to 2016, the state government processed the coal company reorganization and resource product optimizing program [126]. This program aimed to optimize allocation of resources and to improve safety systems. Most of the illegal mining sites have been shut down. Since then, it has begun to develop in a stable stage (Figure 14).

This paper also discussed the bare soil, farmland, and grass land (Figure 4), although some researchers [127,128] thought bared soil is debatable in this kind of research. The Mount Wutai area is a single-season crop area and the seasonal effect is significant. Thus, the farmland is bared for more than half a year. In addition, the landscape of research area is the Loess Plateau, which is no different between bared soil and bared farmland. Therefore, it is inaccurate to evaluate the change of bared land by total area. However, this paper provides a method based on remote sensing and spatial-temporal analysis and it focuses on the spatial variation at the pixel level. In other words, despite that measuring bare soil in such areas might be inaccurate, it works at the pixel level. Therefore, it seems increased from 1989 to 1992, but, if put it in the spatial level, the variation area is grassland or bared soil and, according the tables of GDP and population, it could not be influenced by human activities. It shows a seasonal effect. From this point of view, some grassland, farmland (bared farmland, terraced fields), and bare soil frequently change each other in the study area. Therefore, it is inaccurate to describe their change. Overall, bared soil and farmland together show an increasing tendency. The grassland tends to degrade. Water in Mount Wutai keeps stable.

4.2. Comprehensive Heritage Threats Assessment Index

Disturbances in specific areas are complex and interact drastically. These disturbances sometimes show synergistic or, on the contrary, antagonistic effects to the environment. For disturbances themselves, some are positive to the environment, whereas others are disturbances. Focusing on Mount Wutai, the changes and influence of those four main disturbances can be identified directly; however, it is unknown how it works if they are put together. The comprehensive heritage threats assessment index can explain it well, despite some drawbacks. Based on comprehensive heritage threats assessment index mentioned in the methodology, 12 phases of CHATI, from 1987 to 2018, have been generated (Figure 15). In this map, the deep red color stands for extreme high threats and the deep green stands for extreme safety.

To identity the accuracy of CHATI, the research team selected the most recent CHATI map and did the field validation. In the CHATI map of 2018, most of the threatened areas are located in the north part, the northwest part, and the southwest part of Mount Wutai. Through the fieldwork, it was found that the highly threatened area in the north part of Mount Wutai is the hematite and magnetite producing region, so there are huge human activities at the foot of the mountain. For instance, the mountain is blasted to take minerals and the crude ores are crushed and levitated. During the whole process, a large number of tailings are produced. At the same time, although they are in a legal mining area, due to the poor management in that area, there are visible tailing ponds scattered all over the mountain. There is a long history of hematite mines, as well, in the northeast area named Shanyangping. Focusing on the southwest of Mount Wutai, there are two open basins, located in Doucun Town and Rucun Town, which are mainly used for agricultural production. Theoretically, farmland should not contribute high threats to the environment, however, as this paper mentioned above, farmland in this area shows a significant seasonal effect (see this area in different year maps based on Figure 15) and farmland is bared more than half the year. This pushes a high land stress to the environment. Wutai County, the largest built-up area in the study area, is also located there. There are two open-pit mines around Wutai County. Although it has been closed since 2016, it had been continuously expanding for the past three decades. This area also contributes a large number of threat disturbances to the entire study area. Figure 16 shows the CHATI map combined with fieldworks images.

To illustrate the spatial-temporal variation of the CHATI during the last 3 decades, based on the linear regression and an F-test, mentioned in the methodology, the θ_slope and F maps have been derived (Figure 17). If θ_slope is positive, it indicates that the research object changes with time, and the larger the value, the more obvious the upward trend; otherwise, the research object declines. The interannual rate of change (slope value) of Mount Wutai, from 1987 to 2018, ranges from −0.047 to 0.057. The obvious changes of the research area are located in the third part. Firstly, in Wutai Town and the surrounding villages, there is a significance change in θ_slope map and ‘significance decrease’ is displayed in the F map, because the rapid urbanization and economic development have driven a large number of built-up areas. Secondly, the bared soil and farmland in the southwest display positive in the θ_slope map and ‘increase’ in the F map. After reviewing each selected year’s LULC maps and investigating the local residents around there, there were no significance changes during those years. However, a lot of farmland change between abandonment and redevelopment, combined with the seasonal effect, led to the significant change this area. That is one of the drawbacks of this research. Thirdly, concentrating on the northeast of map which is the core protection area of Mount Wutai, it displays positive in the θ_slope map and shows ‘significance increase’ in the F map. According to the documents [119], from 2000 to 2009 Mount Wutai development rapidly. After that, the local government lead a resettlement plan to protect the core culture landscape. Many houses were demolished and serval settlement village were rebuilt around the core protect area. This results in large and continuous changes. Overall, the comprehensive risk of Wutai Mountain is mostly stable.

5. Conclusions

This research used multi-source (Landsat, MODIS vegetation products, and ASTER-DEM) time-series remote sensing data to identify the threats of the Mount Wutai world landscape heritage area. Firstly, LULCC and human activity disturbances of Mount Wutai were addressed. It shows the forest coverage of Mount Wutai has declined compared to 30 years ago, despite that the local government has implemented several programs to protect the forest and to plant trees, but it still cannot fulfil the deficit of the forest which was destroyed in the middle term. However, the forest coverage has increased in recent years. The most significant change is the built-up area, especially from 1998 to 2004, where the total area of built-up area was doubled. At the same time, due to the high demands of the energy market and rapid economic development, the total mining area has increased fivefold. These dramatic increases directly led to high land stress at the regional level. Secondly, this research proposed a Comprehensive Heritage Area Threats Index (CHATI). This method uses four key normalized indexes constructing a comprehensive system to assess the threats of entire WHA. The accuracy test, based on field validation, has proved this framework works, but there are some drawbacks. Based on the result, five high threat areas have been investigated, and find out the area is high consistent with the LULCC. From this point of view, the CHATI proposed an efficient way to measure the threats of WHA. Thirdly, linear regression and an F-test were adopted in this research to analyze the remarkable characteristics of spatial temporal variation. From the calculation results, the core area of Mount Wutai changed dramatically and trended to less threats. On the contrary, the area related to built-up area and mining showed a high threat.

All these results and analysis has proved reliable and efficient, however, there are some drawbacks need to be discussed here. The time-series Landsat data only use one phase for each selected year. Although this research tried to incorporate the same seasonal data, their seasonal effects were not reflected by the result of LULCC and CHATI. At the same time, 30 years and 12 phases could model the trend of variation, but omits some short-term changes. Therefore, dense time-series datasets can solve this problem, but may also lead to a massive amount of computation. Another shortcoming is with respect to the study area. At specific times, bared soil and farmland (abandoned farmland, terraced fields) have the same features on remote sensing images and it is hard to differentiate accurately, even with DEM, slope, and aspect features. This problem can be solved by incorporating free and open remote sensing datasets to extract useful features to better characterize the LULCC. The last disadvantage is that the pollution indicator has been ignored for CHATI. In the future, we would include PM2.5, rainfall data, etc., to further increase the accuracy and reduce the potential uncertainty.

In conclusion, the study shows that there is considerable potential to apply remote sensing and GIS techniques, in combination with an assessment system, to identify the threat of the World Heritage Area. Long-term time series satellite sensor data could provide rich information to support WHA monitoring and management. The threat and disturbances that provide a new way to identify and monitor the component of WHA were explored. Indeed, there is still the potential to develop this assessment framework, such as the classifier method, the data used, and the indicators that have been ignored. Hopefully this research can raise the public’s awareness on protecting the WHA and may provide an efficient tool to guide public policy making. Finally, we really appreciate the free and open datasets (Landsat, Modis, Modis vegetation product and ASTER-DEM) which enable us to accomplish such detailed research. Without those datasets, there would be no previous research to study. In addition, the whole research area of Mount Wutai cannot be captured for each selected year.

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

[Image omitted. See PDF.]

	Jan	Feb	Mar	Apr	May	June	July	Aug	Sept	Oct	Nov	Dec
2018											OLI
2016									OLI
2013						OLI
2010									TM
2007									TM
2004						TM
2001							TM
1998											TM
1995					TM
1992							TM
1989									TM
1987									TM

Index	Equation	Variables Explain	References
NDVI	NDVI=(NIR−R)/(NIR+R)	NIR: near-infrared band R: infrared band MIR: mid-infrared band BLUE: blue band GREEN: green band SWIR: short-wave infrared band ℓ : soil adjust factor, data range is 0~1, 0 mean extremely high vegetation coverage, contrast 1 is very low. Normally the data is 0.5	[92,93]
MNDWI	MNDWI =Green−MIRGreen+MIR		[46,94]
NSDI	NDSI = (SI + IBI)/2		[95,96,97]
SI	SI=(MIR+R)−(NIR+BULE)(MIR+R)+(NIR+BULE)
IBI	IBI=NDBI−(SAVI+MNDWI)/2NDBI+(SAVI+MNDWI)/2		[98]
SAVI	SI=(SWIR−R)(1+l)SWIR+R+l		[99]

Indicators	Biological Richness Index	Vegetation Coverage Index	Water Network Denseness Index	Land Stress Index	Pollution Load Index	Environmental Restriction Index
Weight	0.35	0.25	0.15	0.15	0.10	Obligatory Target

Parameter	Weight	Indicator
		Indicator	Weight	LULC	Weight
Biological Richness Index	0.35	BI	0
		HQ	1	Forest	0.35
				Grass	0.21
				Water	0.28
				Farmland	0.11
				Built-up Area	0.04
				Unused Land	0.01
Vegetation Coverage Index	0.35	NDVI
Water Network Denseness Index	0.15	NDWI
Land Stress Index	0.15	NDSI

Year		1987	1989	1992	1995	1998	2001	2004	2007	2010	2013	2016	2018
Accuracy	Overall Accuracy (%)	90.47	89.99	93.10	86.50	92.90	90.88	95.17	92.83	93.54	93.16	93.89	94.44
	Kappa Coefficient	0.91	0.91	0.92	0.85	0.93	0.90	0.95	0.92	0.93	0.92	0.93	0.94

[Image omitted. See PDF.]

[Image omitted. See PDF.]

	2011		2012		2013		2014
	The tax amount	Percentage of all tax	The tax amount	Percentage of all tax	The tax amount	Percentage of all tax	The tax amount	Percentage of all tax
Magnesium industry	1438	5.6%	2007	7.8%	2714	8.6%	1551	4.6%
Coal industry	No data	No data	No data	No data	434	1.4%	5428	15.9%
Electric power industry	4060	25.5%	8080	31.3%	8359	26.5%	1958	7.2%

Author Contributions

Conceptualization, X.B., P.Z., and P.T.; methodology, X.B.; validation, X.B., S.G., and P.D.; formal analysis, X.B.; investigation, S.G., C.L., and P.Z.; data curation, S.G.; writing—original draft preparation, X.B., P.D.; writing—review and editing, S.G., P.D., X.B., P.Z., P.T., C.L., and C.Z.; visualization, S.G.; funding acquisition, P.D.

Funding

This study is supported by Open Research Fund of Key Laboratory of Digital Earth Science, Institute of Remote Sensing and Digital Earth, Chinese Academy of Sciences (GrantNo. 2015LDE011), the National Nature Science Foundation of China (GrantNo. 41631176). China Scholarship Council Fund and Newton Fund of British Council (File No. 201806190301).

Acknowledgments

The authors would like to thank Erzhu Li, Jieqiong Luo, Jike Chen, Hongrui Zheng, and Yaping Meng for their help in data processing, and Chuanlong Lin, Ruochen Shang, Zhe Hao, Yun Li, and Yaoyao Wang for their help.

Conflicts of Interest

The authors declare no conflict of interest.

© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

References

1. Poria, Y.; Reichel, A.; Cohen, R. Tourists perceptions of World Heritage Site and its designation. Tour. Manag. 2013, 35, 272-274.

2. Drost, A. Developing sustainable tourism for world heritage sites. Ann. Tour. Res. 1996, 23, 479-484.

3. Leask, A.; Fyall, A. Managing World Heritage Sites; Routledge: Abingdon, UK, 2006.

4. Rössler, M. World heritage cultural landscapes: A UNESCO flagship programme 1992-2006. Landsc. Res. 2006, 31, 333-353.

5. UNESCO World Heritage Centre. The Criteria for Selection; UNESCO World Heritage Centre: Paris, France, 2011.

6. Agapiou, A.; Lysandrou, V.; Alexakis, D.D.; Themistocleous, K.; Cuca, B.; Argyriou, A.; Sarris, A.; Hadjimitsis, D.G. Cultural heritage management and monitoring using remote sensing data and GIS: The case study of Paphos area, Cyprus. Comput. Environ. Urban Syst. 2015, 54, 230-239.

7. Cowley, D.C. Remote Sensing for Archaeological Heritage Management; Europae Archaeologia Consilium: Brussels, Belgium, 2011.

8. Watson, J.E.M.; Dudley, N.; Segan, D.B.; Hockings, M. The performance and potential of protected areas. Nature. 2014, 515, 67-73.

9. Voogt, J.A.; Oke, T.R. Thermal remote sensing of urban climates. Remote Sens. Environ. 2003, 86, 370-384.

10. Pappu, S.; Akhilesh, K.; Ravindranath, S.; Raj, U. Applications of satellite remote sensing for research and heritage management in Indian prehistory. J. Archaeol. Sci. 2010, 37, 2316-2331.

11. Nagendra, H.; Lucas, R.; Honrado, J.P.; Jongman, R.H.G.; Tarantino, C.; Adamo, M.; Mairota, P. Remote sensing for conservation monitoring: Assessing protected areas, habitat extent, habitat condition, species diversity, and threats. Ecol. Indic. 2013, 33, 45-59.

12. Willis, K.S. Remote sensing change detection for ecological monitoring in United States protected areas. Biol. Conserv. 2015, 182, 233-242.

13. Rykiel, E.J. Towards a definition of ecological disturbance. Aust. J. Ecol. 1985, 10, 361-365.

14. Brimblecombe, P.; Grossi, C.M.; Harris, I. Climate change critical to cultural heritage. In Survival and Sustainability; Springer: Berlin, Germany, 2010; pp. 195-205.

15. West, P.; Igoe, J.; Brockington, D. Parks and peoples: The social impact of protected areas. Annu. Rev. Anthropol. 2006, 35, 251-277.

16. Frey, B.S.; Steiner, L. World Heritage List: Does it make sense? Int. J. Cult. Policy 2011, 17, 555-573.

17. Gross, J.E.; Goetz, S.J.; Cihlar, J. Application of remote sensing to parks and protected area monitoring: Introduction to the special issue. Remote Sens. Environ. 2009, 113, 1343-1345.

18. Ali, S.H. Mining, the Environment, and Indigenous Development Conflicts; University of Arizona Press: Tucson, AZ, USA, 2009.

19. Rice, J.; Trujillo, V.; Jennings, S.; Hylland, K.; Hagstrom, O.; Astudillo, A.; Jensen, J.N. Guidance on the Application of the Ecosystem Approach to Management of Human Activities in the European Marine Environment. Available online: https://www.researchgate.net/publication/230659987_Guidance_on_the_Application_of_the_Ecosystem_Approach_to_Management_of_Human_Activities_in_the_European_Marine_Environment (accessed on 30 May 2019).

20. Xu, X.; Yang, G.; Tan, Y.; Zhuang, Q.; Li, H.; Wan, R.; Su, W.; Zhang, J. Ecological risk assessment of ecosystem services in the Taihu Lake Basin of China from 1985 to 2020. Sci. Total Environ. 2016, 554, 7-16.

21. Maekawa, M.; Lanjouw, A.; Rutagarama, E.; Sharp, D. Mountain gorilla tourism generating wealth and peace in post-conflict R wanda. Nat. Resour. Forum 2013, 37, 127-137.

22. Kennedy, R.; Yang, Z.; Braaten, J.; Neldon, P.; Cohen, W. Monitoring landscape dynamics of National Parks in the western United States. In Remote Sensing of Protected Lands; CRC Press: Boca Raton, FL, USA, 2012; pp. 57-75.

23. Wiens, J.; Sutter, R.; Anderson, M.; Blanchard, J.; Barnett, A.; Aguilar-Amuchastegui, N.; Avery, C.; Laine, S. Selecting and conserving lands for biodiversity: The role of Remote Sensing. Remote Sens. Environ. 2009, 113, 1370-1381.

24. Dennison, P.E.; Nagler, P.L.; Hultine, K.R.; Glenn, E.P.; Ehleringer, J.R. Remote monitoring of tamarisk defoliation and evapotranspiration following saltcedar leaf beetle attack. Remote Sens. Environ. 2009, 113, 1462-1472.

25. Jones, J.W.; Hall, A.E.; Foster, A.M.; Smith, T.J. Wetland fire scar monitoring and analysis using archival Landsat data for the Everglades. Fire Ecol. 2013, 9, 133-150.

26. Munns, W.J.E. Society. Assessing risks to wildlife populations from multiple stressors: Overview of the problem and research needs. Ecol. Soc. 2006, 11, 23.

27. Feng, Y.; Liu, Y.; Liu, Y. Assessment, R. Spatially explicit assessment of land ecological security with spatial variables and logistic regression modeling in Shanghai, China. Stoch. Environ. Res. Risk Assess. 2017, 31, 2235-2249.

28. DeFries, R.; Hansen, A.; Newton, A.C.; Hansen, M.C. Increasing isolation of protected areas in tropical forests over the past twenty years. Ecol. Appl. 2005, 15, 19-26.

29. Ingram, J.C.; Dawson, T.P.; Whittaker, R.J. Mapping tropical forest structure in southeastern Madagascar using remote sensing and artificial neural networks. Remote Sens. Environ. 2005, 94, 491-507.

30. Nunes, M.C.; Vasconcelos, M.J.; Pereira, J.M.; Dasgupta, N.; Alldredge, R.J.; Rego, F.C. Land cover type and fire in Portugal: Do fires burn land cover selectively? Landsc. Ecol. 2005, 20, 661-673.

31. Fuller, D.O. Tropical forest monitoring and remote sensing: A new era of transparency in forest governance? Singap. J. Trop. Geogr. 2006, 27, 15-29.

32. Nagendra, H.; Pareeth, S.; Sharma, B.; Schweik, C.M.; Adhikari, K.R. Forest fragmentation and regrowth in an institutional mosaic of community, government and private ownership in Nepal. Landsc. Ecol. 2008, 23, 41-54.

33. Nagendra, H.; Rocchini, D.; Ghate, R.; Sharma, B.; Pareeth, S. Assessing plant diversity in a dry tropical forest: Comparing the utility of Landsat and IKONOS satellite images. Remote Sens. 2010, 2, 478-496.

34. Luo, J.; Du, P.; Samat, A.; Xia, J.; Che, M.; Xue, Z. Spatiotemporal Pattern of PM2.5 Concentrations in Mainland China and Analysis of Its Influencing Factors using Geographically Weighted Regression. Sci. Rep. 2017, 7, 40607.

35. Du, P.; Xia, J.; Zhang, W.; Tan, K.; Liu, Y.; Liu, S. Multiple classifier system for remote sensing image classification: A review. Sensors 2012, 12, 4764-4792.

36. Du, P.; Samat, A.; Waske, B.; Liu, S.; Li, Z. Random forest and rotation forest for fully polarized SAR image classification using polarimetric and spatial features. ISPRS J. Photogramm. Remote Sens. 2015, 105, 38-53.

37. Du, P.; Liu, S.; Xia, J.; Zhao, Y. Information fusion techniques for change detection from multi-temporal remote sensing images. Inf. Fusion 2013, 14, 19-27.

38. Pôças, I.; Cunha, M.; Pereira, L.S. Remote sensing based indicators of changes in a mountain rural landscape of Northeast Portugal. Appl. Geogr. 2011, 31, 871-880.

39. Gitelson, A.A.; Merzlyak, M.N. Signature analysis of leaf reflectance spectra: Algorithm development for remote sensing of chlorophyll. J. Plant Physiol. 1996, 148, 494-500.

40. Liu, Y.; Sarnat, J.A.; Kilaru, V.; Jacob, D.J.; Koutrakis, P. Estimating ground-level PM2. 5 in the eastern United States using satellite Remote Sensing. Environ. Sci. Technol. 2005, 39, 3269-3278.

41. Van Donkelaar, A.; Martin, R.V.; Park, R.J. Estimating ground-level PM2. 5 using aerosol optical depth determined from satellite Remote Sens. J. Geophys. Res. Atmos. 2006, 111.

42. Kimes, D.; Sellers, P.J. Inferring hemispherical reflectance of the Earth's surface for global energy budgets from remotely sensed nadir or directional radiance values. Remote Sens. Environ. 1985, 18, 205-223.

43. Bastiaanssen, W.G.; Menenti, M.; Feddes, R.; Holtslag, A.A.M. A remote sensing surface energy balance algorithm for land (SEBAL). 1. Formulation. J. Hydrol. 1998, 212, 198-212.

44. Jackson, T.J., III. Measuring surface soil moisture using passive microwave Remote Sensing. Hydrol. Process. 1993, 7, 139-152.

45. Weng, Q. A remote sensing? GIS evaluation of urban expansion and its impact on surface temperature in the Zhujiang Delta, China. Int. J. Remote Sens. 2001, 22, 1999-2014.

46. Gao, B.-C. NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sens. Environ. 1996, 58, 257-266.

47. Bai, Z.; Dent, D.; Olsson, L.; Schaepman, M. Global Assessment of Land Degradation and impRovement: 1. Identification by Remote Sensing; ISRIC-World Soil Information: Wageningen, The Netherlands, 2008.

48. Curran, P. Multispectral remote sensing of vegetation amount. Prog. Phys. Geogr. 1980, 4, 315-341.

49. Lu, D. The potential and challenge of remote sensing-based biomass estimation. Int. J. Remote Sens. 2006, 27, 1297-1328.

50. Ruimy, A.; Saugier, B.; Dedieu, G. Methodology for the estimation of terrestrial net primary production from remotely sensed data. J. Geophys. Res. Atmos. 1994, 99, 5263-5283.

51. Field, C.B.; Randerson, J.T.; Malmström, C.M. Global net primary production: Combining ecology and Remote Sensing. Remote Sens. Environ. 1995, 51, 74-88.

52. Lefsky, M.A.; Cohen, W.B.; Parker, G.G.; Harding, D. J Lidar remote sensing for ecosystem studies: Lidar, an emerging remote sensing technology that directly measures the three-dimensional distribution of plant canopies, can accurately estimate vegetation structural attributes and should be of particular interest to forest, landscape, and global ecologists. BioScience 2002, 52, 19-30.

53. Dandois, J.P.; Ellis, E.C. Remote sensing of vegetation structure using computer vision. Remote Sens. 2010, 2, 1157-1176.

54. Kerr, J.T.; Ostrovsky, M. From space to species: Ecological applications for Remote Sensing. Trends Ecol. Evol. 2003, 18, 299-305.

55. Quattrochi, D.A.; Luvall, J.C. Thermal infrared remote sensing for analysis of landscape ecological processes: Methods and applications. Landsc. Ecol. 1999, 14, 577-598.

56. Van Donkelaar, A.; Martin, R.V.; Brauer, M.; Hsu, N.C.; Kahn, R.A.; Levy, R.C.; Lyapustin, A.; Sayer, A.M.; Winker, D.M. Global estimates of fine particulate matter using a combined geophysical-statistical method with information from satellites, models, and monitors. Environ. Sci. Technol. 2016, 50, 3762-3772.

57. Zheng, Y.; Zhang, Q.; Liu, Y.; Geng, G.; He, K. Estimating ground-level PM2. 5 concentrations over three megalopolises in China using satellite-derived aerosol optical depth measurements. Atmos. Environ. 2016, 124, 232-242.

58. Guo, Y.; Tang, Q.; Gong, D.-Y.; Zhang, Z. Estimating ground-level PM2. 5 concentrations in Beijing using a satellite-based geographically and temporally weighted regression model. Remote Sens. Environ. 2017, 198, 140-149.

59. Simic, A.; Chen, J.M.; Liu, J.; Csillag, F. Spatial scaling of net primary productivity using subpixel information. Remote Sens. Environ. 2004, 93, 246-258.

60. Bian, J.-H.; Li, A.; Song, M.; Ma, L.; Jiang, J. Reconstruction of NDVI time-series datasets of MODIS based on Savitzky-Golay filter. J. Remote Sens. 2010, 14, 725-741.

61. Lu, D.; Xu, X.; Tian, H.; Moran, E.; Zhao, M.; Running, S. The effects of urbanization on net primary productivity in southeastern China. Environ. Manag. 2010, 46, 404-410.

62. Xu, H.; Wang, M.; Shi, T.; Guan, H.; Fang, C.; Lin, Z. Prediction of ecological effects of potential population and impervious surface increases using a remote sensing based ecological index (RSEI). Ecol. Indic. 2018, 93, 730-740.

63. Srinivasa Rao, Y.; Jugran, D.K. Delineation of groundwater potential zones and zones of groundwater quality suitable for domestic purposes using remote sensing and GIS. Hydrol. Sci. J. 2003, 48, 821-833.

64. Miller, R.B.; Small, C. Cities from space: Potential applications of remote sensing in urban environmental research and policy. Environ. Sci. Policy 2003, 6, 129-137.

65. Luo, J.; Zhou, T.; Du, P.; Xu, Z. Spatial-temporal variations of natural suitability of human settlement environment in the Three Gorges Reservoir Area—A case study in Fengjie County, China. Front. Earth Sci. 2019, 13, 1-17.

66. Yinkang, C.F.; Feng, Z. Analysis of the spatial distribution pattern of urban land price with geostatistics. J. Naijing Univ. Nat. Sci. 1999, 6, 719-723.

67. Luo, J.; Du, P.; Samat, A.; Feng, L. Evaluation on the natural suitability of urban human settlement environment using multisource data. In Proceedings of the 2015 Joint Urban Remote Sensing Event (JURSE), Lausanne, Switzerland, 30 March-1 April 2015; pp. 1-4.

68. Patino, J.E.; Duque, J.C.; Pardo-Pascual, J.E.; Ruiz, L.A. Using remote sensing to assess the relationship between crime and the urban layout. Appl. Geogr. 2014, 55, 48-60.

69. Xu, Z.; Li, Q. Integrating the empirical models of benchmark land price and GIS technology for sustainability analysis of urban residential development. Habitat Int. 2014, 44, 79-92.

70. Xu, Y.; Sun, J.; Zhang, J.; Xu, Y.; Zhang, M.; Liao, X. Combining AHP with GIS in synthetic evaluation of environmental suitability for living in China's 35 major cities. Int. J. Geogr. Inf. Sci. 2012, 26, 1603-1623.

71. Ouzounis, G.K.; Syrris, V.; Pesaresi, M. Multiscale quality assessment of Global Human Settlement Layer scenes against reference data using statistical learning. Pattern Recognit. Lett. 2013, 34, 1636-1647.

72. Shen, L.; Kyllo, J.; Guo, X. An integrated model based on a hierarchical indices system for monitoring and evaluating urban sustainability. Sustainability 2013, 5, 524-559.

73. Xu, C.; Liu, M.; An, S.; Chen, J.; Yan, P. Assessing the impact of urbanization on regional net primary productivity in Jiangyin County, China. J. Environ. Manag. 2007, 85, 597-606.

74. Wu, C.; Murray, A.T. Estimating impervious surface distribution by spectral mixture analysis. Remote Sens. Environ. 2003, 84, 493-505.

75. Du, H.; Cai, W.; Xu, Y.; Wang, Z.; Wang, Y.; Cai, Y. Quantifying the cool island effects of urban green spaces using remote sensing Data. Urban For. Urban Green. 2017, 27, 24-31.

76. Du, P.; Bai, X.; Luo, J.; Li, E.; Lin, C. Advances of urban Remote Sensing. Nanjing Xinxi Gongcheng Daxue Xuebao 2018, 10, 16-29.

77. Li, S.; Zhao, G.; Wilde, S.A.; Zhang, J.; Sun, M.; Zhang, G.; Dai, L. Deformation history of the Hengshan-Wutai-Fuping Complexes: Implications for the evolution of the Trans-North China Orogen. Gondwana Res. 2010, 18, 611-631.

78. Liu, S.; Pan, Y.; Xie, Q.; Zhang, J.; Li, Q. Archean geodynamics in the Central Zone, North China Craton: Constraints from geochemistry of two contrasting series of granitoids in the Fuping and Wutai complexes. Precambrian Res. 2004, 130, 229-249.

79. Stevenson, D. Visions of Mañjuśrī on Mount Wutai. Relig. China Pract. 1996, 3, 203.

80. Lin, W.-C. Building a Sacred Mountain: The Buddhist Architecture of China's Mount Wutai; University of Washington Press: Seattle, WA, USA, 2014.

81. Ryan, C.; Gu, H. Constructionism and culture in research: Understandings of the fourth Buddhist Festival, Wutaishan, China. Tour. Manag. 2010, 31, 167-178.

82. Wang, K.; Li, J.; Hao, J.; Li, J.; Shaoping, Z. The Wutaishan orogenic belt within the Shanxi Province, northern China: A record of late Archaean collision tectonics. Precambrian Res. 1996, 78, 95-103.

83. Wutai Bureau of Statistics of Wutai County. Yearbook; Wutai Bureau of Statistics of Wutai County: Xinzhou, China, 2010.

84. Wutai Bureau of Statistics of Wutai County. Yearbook; Wutai Bureau of Statistics of Wutai County: Xinzhou, China, 2013.

85. Zhang, W.; Wang, J. School of Tourism Management, Xiangtan University; On Development of Red Tourism Resources in China—Taking Wutai County in Shanxi Province as an Example. J. Hunan Univ. Technol. Soc. Sci. Ed. 2013, 4, 9.

86. Thenkabail, P.S.; Schull, M.; Turral, H. Ganges and Indus river basin land use/land cover (LULC) and irrigated area mapping using continuous streams of MODIS data. Remote Sens. Environ. 2005, 95, 317-341.

87. Yang, X.; Huang, Y.; Dong, P.; Jiang, D.; Liu, H. An updating system for the gridded population database of China based on remote sensing, GIS and spatial database technologies. Sensors 2009, 9, 1128-1140.

88. Mengistu, D.A.; Salami, A.T. Technology. Application of remote sensing and GIS inland use/land cover mapping and change detection in a part of south western Nigeria. Afr. J. Environ. Sci. Technol. 2007, 1, 99-109.

89. Alo, C.A.; Pontius, R.G., Jr. Identifying systematic land-cover transitions using remote sensing and GIS: The fate of forests inside and outside protected areas of Southwestern Ghana. Environ. Plan. B Plan. Des. 2008, 35, 280-295.

90. Liu, J.; Linderman, M.; Ouyang, Z.; An, L.; Yang, J.; Zhang, H. Ecological degradation in protected areas: The case of Wolong Nature Reserve for giant pandas. Science 2001, 292, 98-101.

91. Kairo, J.G.; Kivyatu, B.; Koedam, N. Application of Remote Sensing and GIS in the Management of Mangrove Forests Within and Adjacent to Kiunga Marine Protected Area, Lamu, Kenya. Environ. Dev. Sustain. 2002, 4, 153-166.

92. Running, S.W. Estimating terrestrial primary productivity by combining remote sensing and ecosystem simulation. In Remote Sensing of Biosphere Functioning; Springer: Berlin, Germany, 1990; pp. 65-86.

93. Myneni, R.B.; Hall, F.G.; Sellers, P.J.; Marshak, A.L. The interpretation of spectral vegetation indexes. IEEE Trans. Geosci. Remote Sens. 1995, 33, 481-486.

94. Han-Qiu, X.U. A Study on Information Extraction of Water Body with the Modified Normalized Difference Water Index (MNDWI). J. Remote Sens. 2005, 5, 589-595.

95. Zhou, C.-L.; Xu, H.-Q. A Spectral Mixture Analysis and Mapping of Impervious Surfaces in Built-up Land of Fuzhou City. J. Image Graph. 2007, 5, 875-881.

96. Xu, H. Assessment of ecological change in soil loss area using remote sensing technology. Trans. Chin. Soc. Agric. Eng. 2013, 29, 91-97.

97. Xu, H. Fast information extraction of urban built-up land based on the analysis of spectral signature and normalized difference index. Geogr. Res. 2005, 24, 311-320.

98. Xu, H. A new index for delineating built-up land features in satellite imagery. Int. J. Remote Sens. 2008, 29, 4269-4276.

99. Huete, A.R. A soil-adjusted vegetation index (SAVI). Remote Sens. Environ. 1988, 25, 295-309.

100. Zhu, G.; Liu, Y.; Ju, W.; Chen, J. Evaluation of topographic effects on four commonly used vegetation indices. Yaogan Xuebao-J. Remote Sens. 2013, 17, 210-234.

101. Ministry of Ecology and Environment of the P.R.China. Technical Specifications for Evaluation of Ecological Environment Conditions; Ministry of Ecology and Environment of the P.R.China, Ed.; 2015.

102. Seber, G.A.; Lee, A.J. Linear Regression Analysis; John Wiley & Sons: Hoboken, NJ, USA, 2012; Volume 329.

103. Lix, L.M.; Keselman, J.C.; Keselman, H. Consequences of assumption violations revisited: A quantitative review of alternatives to the one-way analysis of variance F test. Rev. Educ. Res. 1996, 66, 579-619.

104. Tripathy, G.; Ghosh, T.; Shah, S.D. Monitoring of desertification process in Karnataka state of India using multi-temporal remote sensing and ancillary information using GIS. Int. J. Remote Sens. 1996, 17, 2243-2257.

105. Seto, K.C.; Kaufmann, R.K. Modeling the drivers of urban land use change in the Pearl River Delta, China: Integrating remote sensing with socioeconomic data. Land Econ. 2003, 79, 106-121.

106. Johnson, L.; Bosch, D.; Williams, D.; Lobitz, B.M. Remote sensing of vineyard management zones: Implications for wine quality. Appl. Eng. Agric. 2001, 17, 557.

107. Johnson, B.A.; Iizuka, K. Integrating OpenStreetMap crowdsourced data and Landsat time-series imagery for rapid land use/land cover (LULC) mapping: Case study of the Laguna de Bay area of the Philippines. Appl. Geogr. 2016, 67, 140-149.

108. Belgiu, M.; Drǎguţ, L. Comparing supervised and unsupervised multiresolution segmentation approaches for extracting buildings from very high resolution imagery. ISPRS J. Photogramm. Remote Sens. 2014, 96, 67-75.

109. Al-Saady, Y.; Merkel, B.; Al-Tawash, B.; Al-Suhail, Q. Land use and land cover (LULC) mapping and change detection in the Little Zab River Basin (LZRB), Kurdistan Region, NE Iraq and NW Iran. FOG-Freib. Online Geosci. 2015, 43, 1-32.

110. Cohen, J. A coefficient of agreement for nominal scales. Educ. Psychol. Meas. 1960, 20, 37-46.

111. Thomlinson, J.R.; Bolstad, P.V.; Cohen, W.B. Coordinating methodologies for scaling landcover classifications from site-specific to global: Steps toward validating global map products. Remote Sens. Environ. 1999, 70, 16-28.

112. Zhang, M. China's Poor Regions: Rural-Urban. Migration, Poverty, Economic Reform and Urbanisation; Routledge: Abingdon, UK, 2004.

113. Shanxi Bureau of statistics of Shanxi Province. Yearbook; Shanxi Bureau of statistics of Shanxi Province: Shanxi, China, 2018.

114. Wutai Bureau of statistics of Wutai County. Yearbook; Wutai Bureau of Statistics of Wutai County: Xinzhou, China, 2005.

115. Myneni, R.; Tucker, C.; Asrar, G.; Keeling, C.D. Interannual variations in satellite-sensed vegetation index data from 1981 to 1991. J. Geophys. Res. Atmos. 1998, 103, 6145-6160.

116. Wang, Y.P.; Murie, A. The process of commercialisation of urban housing in China. Urban Stud. 1996, 33, 971-989.

117. Chen, J.; Guo, F.; Wu, Y. One decade of urban housing reform in China: Urban housing price dynamics and the role of migration and urbanization, 1995-2005. Habitat Int. 2011, 35, 1-8.

118. Logan, J.R.; Fang, Y.; Zhang, Z. The winners in China's urban housing reform. Hous. Stud. 2010, 25, 101-117.

119. Wutai Bureau of statistics of Wutai County. Yearbook; Wutai Bureau of Statistics of Wutai County: Xinzhou, China, 2018.

120. Wutai Bureau of statistics of Wutai County. Yearbook; Wutai Bureau of Statistics of Wutai County: Xinzhou, China, 2014.

121. Xue, L.; Wang, M.Y.; Xue, T. ‘Voluntary'poverty alleviation resettlement in China. Dev. Chang. 2013, 44, 1159-1180.

122. Tu, J. Coal mining safety: China's Achilles' heel. China Secur. 2007, 3, 36-53.

123. Cao, X. Policy and regulatory responses to coalmine closure and coal resources consolidation for sustainability in Shanxi, China. J. Clean. Prod. 2017, 145, 199-208.

124. Gao, H.Q.; Xiong, Y.Z.; Hu, S.W. World Financial Crisis and Its Impact on China. Stud. Int. Financ. 2008, 11.

125. Zhang, J.; Fu, M.; Geng, Y.; Tao, J. Energy saving and emission reduction: A project of coal-resource integration in Shanxi Province, China. Energy Policy 2011, 39, 3029-3032.

126. Xie, H.; Wang, J.; Shen, B.; Liu, J.; Jiang, P.; Zhou, H.; Liu, H.; Wu, G. New idea of coal mining: Scientific mining and sustainable mining capacity. J. China Coal Soc. 2012, 37, 1069-1079.

127. Metternicht, G.; Zinck, J.A. Remote sensing of soil salinity: Potentials and constraints. Remote Sens. Environ. 2003, 85, 1-20.

128. Fang, B.; Lakshmi, V. Soil moisture at watershed scale: Remote sensing techniques. J. Hydrol. 2014, 516, 258-272.