Ecosystem Health Assessment of World Natural Heritage Sites Based on Remote Sensing and Field Sampling Veriﬁcation: Bayanbulak as Case Study

: Monitoring the ecosystem health for world natural heritage sites is essential for protecting them and beneﬁts the formulation of more targeted protection policies. This study used Bayanbulak world natural heritage site as a case, established a framework for assessing the ecosystem health through remote sensing based on the parameters of ecosystem vigour, organization, resilience, and services. Then, we veriﬁed the obtained results through ﬁeld sampling. The results show that the ecosystem health in the overall study area had declined over time, however, the health within the property zone remained at high levels and stable. The area proportion of low health was low and primarily distributed in the bu ﬀ er zone. Thus, in general, the ecosystem in the study area was healthy. Besides, the ecosystem health exhibited distinct spatial agglomeration characteristics, and the degree of agglomeration enhanced over time. In addition, the ﬁeld vegetation samplings were consistent with the changes in the ecosystem health levels, therefore, the result of RS monitoring of ecosystem health were credible. Thus, this study provides a scientiﬁc basis for heritage managers to formulate suitable ecological protection policies and should aid further research on the ecological monitoring of heritage sites.


Introduction
World natural heritage (WNH) sites are natural areas with global outstanding universal value, and serve as regional and background references for mankind's cognition with respect to evolution of natural, biodiversity, bioecological processes, and natural landscape beauty [1]. Thus, they are of great significance for protecting biodiversity, maintaining ecological health, and striking a balance between human development and the environment [2]. However, WNH sites face many threats from human activities and natural disasters, such as development and construction, the overutilization of resources, environmental pollution, pests and diseases, earthquakes, and mudslides [3,4]. Therefore, United Nations Educational, Scientific and Cultural Organization has strengthened the monitoring of WNH sites and added those with significant problems to the "Endangered World Heritage List" to urge the countries where these sites are located to take necessary measures to restore and protect them. By 2018, this list contained 16 WNH sites from around the world [5]. Given these facts, the ecological health assessment and monitoring of WNH sites is attracting a lot of research interest, in order to develop scientific methods for providing early warnings to protect these sites [6,7].

Study Area
Bayanbulak is located in the bottom swampland of the Youerdusi Basin in the middle of Xinjiang Tianshan (Figure 1). It belongs to the Kaidu River basin, and is the upriver catchment area of the river. It is the best representative of the large intermontane basins of the Tianshan Mountains, and is a typical example of an alpine wetland ecosystem in the arid temperate zone. It is also an excellent example of the beautiful landscape of the Tianshan Mountains, consisting of bending rivers and marshes. The topography of Bayanbulak is gentle, with that in the west being slightly higher. The highest altitude is 2600 m, while the lowest is 2390 m, and the vertical relief is only 210 m. It is surrounded by mountains and has a temperate continental arid climate. The annual average temperature is −4.6 • C, and the extreme high and low temperatures are 28.3 • C and −48.1 • C, respectively. The annual average precipitation is 276 mm and mainly occurs from June to August, which account for 50-70% of the total. Snowfall is concentrated between January and March, and the annual average snowfall is 70.5 mm. The annual evaporation is 1128 mm, while the average relative humidity is 69% [40]. The soil in the region is predominantly alpine basin swamp soil. Because of the frosty climate and frozen soil, there are few large trees in the region, with the primary vegetation being typical alpine swamps and meadows. However, abundant species of birds are present, and it is the largest National Nature Reserve (NNR) for swans in China. The study area of Bayanbulak included the property zone and its surrounding buffer zone (Figure 1). The property zone covers the best-preserved areas of the alpine wetland ecosystem and falls with the boundaries of the Bayanbulak NNR. To better protect the integrity of the heritage site, a protective transition zone (i.e., a buffer zone) has been set up around the property zone. In addition, the property zone is divided into a prohibited zone, a restricted zone, and an exhibition zone based on the distribution characteristics and sensitivity of the protected entities [40]. Among these zones, the prohibited zone is the core area where both the heritage value and degree of vulnerability are very high, and all potentially destructive activities are strictly prohibited. On the other hand, the restricted zone is an area with high bioecological value and relatively high ecological vulnerability. While activities related to ecotourism and science education can be carried out here, the originality and integrity of the area must be maintained. Finally, the exhibition zone is where indigenous people are allowed live. In addition, sightseeing by tourists is also permitted. However, all activities must be carried out while ensuring that the ecological environment in the zone is protected.

Data Sources
The data used in this study included spatial data and field sampling data. Spatial data: World natural heritage boundary of Bayanbulak comes from Xinjiang Tianshan World Natural Heritage Declaration [39]; digital elevation model (DEM) (the resolution was 30 meter) was downloaded from Geospatial Data Cloud (http://www.gscloud.cn/); and Landsat thematic mapper/enhanced thematic mapper RS images of the Bayanbulak WNH site for 2000, 2011, and 2018 (the resolution was 30 meter) were obtained from the United States Geological Survey (http://glovis.usgs.gov/). Field sampling data: The data was collected from the field sampling of vegetation in 2018 and 2019.

Data Processing
We used the software packages ENVI 5.3 and ArcGIS 10.5 to pre-process the RS image data, which were subjected to band synthesis, radiometric corrections, atmospheric corrections, mosaic making, and cutting. To made the spatial data matched perfectly, we used ArcGIS 10.5 to define and project the spatial data layers in the same coordinate system (WGS_1984_UTM_45N), and used the spatial analysis function of ArcGIS 10.5 for spatial overlays analysis of EHA results and spatial display of EHA indicators. To obtained ecosystem landscape types, a supervised classification method was used in combination with the DEM and normalised difference vegetation index (NDVI), and the field data in 2018 and 2019 as auxiliary data for manual visual interpretation. The study area was divided into swamp grassland, high-coverage grassland, medium-coverage grassland, and low-coverage grassland, water, riverbed, construction land, sand, and bare land, and the ecosystem landscape types in 2000, 2011, and 2018 were determined. The landscape pattern index was calculated using Fragstats 4.2. The field sampling data were processed using Microsoft Excel. Correlation analyses were performed using SPSS 25.0. Spatial autocorrelation analysis was calculated using the software Geoda.

Ecosystem Health Assessment Framework
Based on previous studies, we divided the research scales of EHA into three levels scale, e.g., large, medium, and small scale. The large-scale research refers to the study with global or national as evaluation units [24,27]; the medium-scale research was taking provincial or municipal administrative district as evaluation units [25,31]; and the small-scale research was based on the evaluation unit of a specific nature conservation site [41]. Our study is the world natural heritage of Bayanbulak, and belongs to the small-scale ecosystem health assessment. We used the VORS model to establish an EHA framework for the Bayanbulak WNH site based on previous studies [12,42,43]. The framework includes the following parts: Establishing an indicators system for EHA, analysing the space-time evolution of the ecosystem health and the effects of spatial agglomeration, and evaluating the accuracy of obtained results by comparing it with field vegetation sampling data ( Figure 2). Note: Ecosystem health assessment (EHA), ecosystem vigour (EV), ecosystem organization (EO), ecosystem resilience (ER), ecosystem service (ES), landscape heterogeneity (LH), landscape connectivity (LC), connectivity of patches with important ecological functionality (IC). RS data is remote sensing image data; landscape types data comes from remote sensing image.

Ecosystem Health Assessment
• Ecosystem Health Assessment Model The VORS model uses the parameters EV, EO, ER, and ES to determine the level of ecosystem health. Owing to the complexity of evaluating these indicators and the differences in the properties and dimensions of the data involved, the various indicators had to be made comparable. For this, different standardisations methods were employed, based on the role of indicators in the ecosystem [44,45].
Indicators with positive correlation effects, the standardized formula is: Indicators with negative correlation effects, the standardized formula is: Based on the actual conditions of the study area, a 1 km × 1 km grid was selected for monitoring and evaluating the ecosystem health in the study area. Since each factor is indispensable for assessing the ecosystem health, the four factors were multiplied. In addition, because the size of the standardised data after the four multiplications was too small, and given that the differences between the data would be reduced after the multiplication process, the fourth root of the product was taken [12,46]. The formula is: where, EH is ecosystem health value; EV is ecosystem vigour; EO is ecosystem organization; ER is ecosystem resilience; ES is ecosystem service value.

•
Ecosystem Health Assessment Indicators The indicators system of EHA as showed Table 1. EV refers to the energy input and nutrient cycle capacity of the ecosystem, and can be expressed as the metabolism or primary productivity of the ecosystem [10]. On the other hand, the NDVI reflects the space-time dynamics of the vegetation distribution and productivity [49]. Previous studies have shown that the NDVI is highly positively correlated with vegetation productivity [50][51][52]. Therefore, we used the NDVI as an indicator of EV. The formula is: where, NIR is the near infrared band; RED is the infrared band. EO represents the structural stability of the ecosystem, describes the interrelationships among its components, and reflects the complexity of its structure, including its diversity and connectivity [11]. Landscape ecology emphasises the interactions between spatial patterns and ecological processes and is of great significance to the ecosystem, providing an effective way of monitoring and assessing its health [21,53]. It is usually measured using the landscape pattern indices, including the landscape heterogeneity (LH), landscape connectivity (LC), and connectivity of patches with important ecological functionality (IC). LH is maintained through landscape diversity [54] and can be quantified using Shannon's diversity index. LC is quantified through landscape fragmentation. Since wetlands play an important role in determining the ecosystem functionality in Bayanbulak, we used the swamp grassland fragmentation index and the patch cohesion index to quantify the IC.
The weight coefficient model was used to quantify the EO [55]. Ecosystem health is affected by LH and LC, which describe different aspects of the ecosystem structure and thus are equally important [12,56]. Moreover, the weights assigned to IC should not be higher than that for the overall patch connectivity [12]. Thus, the weights assigned to LH, LC, and IC were 0.35, 0.35, and 0.3, respectively. The degree of fragmentation is the core metric that determines the connectivity, and the weight assigned to the fragmentation degree should be higher than that for the patch cohesion index [42]. Thus, the swamp grassland fragmentation degree and patch cohesion index were assigned weights of 0.2 and 0.1, respectively. The formula is: where, LH is landscape heterogeneity; LC is landscape connectivity; IC is connectivity of patches with important ecological functionality; SHDI is Shannon's diversity index; FN is overall landscape fragmentation; FN1 is swamp grassland landscape fragmentation; COHESION is swamp grassland patch cohesion index. ER refers to the ability of the ecosystem to maintain its structure and patterns when subjected to external interference and is reflective of its resistance to external disturbances, as well as its ability to recover [57]. A healthy ecosystem should be resilient enough to be able to withstand various small-scale disturbances, and different types of ecosystems will have different ecosystem resilience coefficients. In this study, the ecosystem resilience coefficient was determined based on previous studies and expert knowledge (Table 2) [12,27,42,56]. We used area-weighted ecosystem resilience coefficients for all the ecosystem types quantified. The formula is: where, A i is area ratio of ecosystem type i; RC i is resilience coefficient of ecosystem type i; n is the number of ecosystem type. Finally, ES refers to the life support products and services that are directly or indirectly available because of the structure, processes, and functions of an ecosystem [58]. A healthy ecosystem should have a high ES supply capacity. In this study, we used the value equivalent method to evaluate the ES value; this method distinguishes between the different types of ES functions and quantifiable standards to calculate the value equivalents for the various service functions in different types of ecosystems and then combines them with its area for the entire ecosystem [59,60]. The construction of the value equivalents of the ecological service functions is the key to determining their values. Based on the work of Xie et al. [47,48] and the actual conditions in the study area, we were able to determine the value equivalent for the ES functions in the study area (Table 3). The formula is: where, ES is the total ecosystem service value; A i is the area of ecosystem type i; VC ij is the j type of ecosystem service value of ecosystem type i.

Spatial Autocorrelation Analysis
Spatial autocorrelation analyses were performed to study the spatial dependence and agglomeration patterns of the ecosystem health for Bayanbulak. This included both global and local spatial autocorrelation analyses. Global spatial autocorrelation reveals the overall trends in the global spatial correlation by calculating the global Moran's I index, whose range is [−1,1]. Index values greater than 0 indicate that the spatial correlation is positive; the larger the value, the more positive the correlation and hence the stronger the spatial agglomeration. On the other hand, index values lower than 0 indicate that the spatial correlation is negative; the smaller the value, the more negative the correlation and the greater the spatial difference. Finally, an index value of 0 means that there is no spatial correlation and that the spatial units are randomly distributed. The local spatial autocorrelation reflects the local characteristics of the spatial correlation and is determined by calculating the local Moran's I index (LISA). Index values greater than 0 mean that the spatial units exhibit high-high or low-low agglomeration, while values smaller than 0 mean that the spatial units exhibit high-low or low-high agglomeration. The formulas are listed as followed [61,62]: where, n is the number of units in the study area; x i is the value of the spatial unit i; x j is the value of the spatial unit j; W ij is the spatial matrix of spatial autocorrelation; x is the mean value.

Verification Based on Field Data
In order to further verify the accuracy and effectiveness of the proposed EHA model, we used the field vegetation sampling data to indirectly evaluate the accuracy of the model with respect to the study area. We collected field vegetation samples in the study area in July 2018 and July 2019, based on the accessibility of the roads. Samples from a total of 47 plots were collected. Next, as per previous studies [63][64][65][66], the EHA results were verified based on the vegetation coverage, diversity, richness, and evenness of the field vegetation. We used Spearman rank correlation coefficient to elucidate the relationship between the EHA results and the biodiversity of the field vegetation, with the aim of determining the accuracy of the EHA results.
We selected the Vegetation Coverage, Simpson Diversity Index, Shannon-Wiener Diversity Index, Margalef Richness Index, Pielou Evenness Index as the field vegetation biodiversity. The formula is [67]: Shannon-Wiener diversity index = − P i lnP i Margalef richness index = (S−1)/lnN (12) Pielou evenness index = − P i lnP i /lnS (13) where, P i = N i /N, N i is the number of vegetation species i; N is the total number of sample species; S is the total number of vegetation in the sample.

Changes in Ecosystem Landscape Types
From Figures 3 and 4, we found that swamp grasslands and high-and medium-coverage grassland are the main landscape types in the study area, accounting for 95%, 95%, and 94% of the total area in 2000, 2011, and 2018, respectively; while the proportion of the construction land was relatively low, accounting for only 0.4%, 0.5%, and 0.4% of the study area. Overall, there was no change in the ecosystem landscape types structure from 2000 to 2018; however, there were some changes in the areas of the different types. Specifically, the area of medium-coverage grassland increased significantly, from 14.5% in 2000 to 22% in 2018. In addition, the area of the high-coverage grassland decreased, going from 38.5% in 2000 to 30.1% in 2018. Further, the area of the water reduced by 2.3% over the same period, while those of the other landscape types remained stable.   Table 4. In general, the degree of landscape fragmentation increased in all three years. However, the extent of fragmentation remained small (0.0061-0.0084). Specifically, the fragmentation index for swamp grassland remained almost constant, while that for high-coverage grassland increased initially and then decreased. Finally, that for medium-coverage grassland fell.

Changes in EHA indicators
The spatial distributions of and changes in the EHA indicators in the overall study area during 2000, 2011, and 2018 were determined ( Figure 5, Table 5   In contrast, EO remained essentially unchanged in the study area, as well as the property zone and buffer zone, during the three investigated years. Its overall value was relatively low, and the mean values for the study area, property zone, and buffer zone were 0.54, 0.63, and 0.45, respectively. Its spatial distribution also remained stable, especially within the property zone, where its value was highest, though there was a slight change in its value in the northern part of the buffer zone. The value of ER remained almost stable in the overall study area. Specifically, within the property zone, its value was generally high; it declined slightly from 2000 (0.9024) to 2011 (0.8796), but remained stable from 2011 and 2018. The spatial structure within property zone remained stable. On the other hand, the buffer zone showed a marked decline from 2000 (0.7802) to 2011 (0.7177), and then a slight decrease in 2011-2018. Further, its spatial distribution suggested that it had declined in the northwestern part in 2011 and 2018.
In general, ES declined in the overall study during 2000-2018 and, in particular, from 2000 (0.6507) to 2011 (0.582). Further, within the property zone, which value of ES was high, had a slow trend of decline (the mean value of 2000, 2011, and 2018 was 0.7928, 0.7544, and 0.7445, respectively). However, in the buffer zone, the ES value remained low, and its rate of decline was significant, especially from 2000 (0.4881) to 2011 (0.3826), while in 2018, it fell to 0.3477. In terms of spatial distribution, it remained essentially unchanged within the property zone, while that in the buffer zone changed significantly. Compared with the case in 2000, the value in the northern and northwestern parts declined significantly during 2011, while in 2018, the value in the southwestern and northwestern parts declined significantly.
In summary, all the EHA indicators showed declining trends, with the decline in ES being the most significant. On comparing the ecosystem landscape types and the EHA indicator values for 2000, 2011, and 2018, it can be inferred that the distribution of swamp grassland and high-, medium-, and low-coverage grassland is the primary factor influencing the values of the EHA indicators in the study area.

Changes in Ecosystem Health
The overall ecosystem health in Bayanbulak also showed a declining trend during the investigated periods (Table 5). For instance, it declined from 0.6802 in 2000 to 0.6251 in 2018, with the decline occurring primarily from 2000 (0.6802) to 2011 (0.6324). While the variations within the property zone were small, there was a slight decline here as well, however, the overall health remained stable. Therefore, the changes in the overall study area were mainly caused by the changes within the buffer zone, with the decline from 2000 (0.6046) to 2011 (0.5284) being especially significant, however, the rate of decline became smaller from 2011 (0.5284) to 2018 (0.5075).
In terms of proportions of area in ecosystem health levels ( Figure 6). The proportions of the area with "good" ecosystem health in 2000, 2011, and 2018 were 47%, 44%, and 44%, respectively, thus it remained stable. However, the proportion of the area with "relatively good" health declined significantly, from 47% in 2000 to 26% in 2018; the rate of decline from 2000 to 2011 was 15%. In addition, the proportion of the area with "ordinary" health showed an increase, going from 6% in 2000 to 27% in 2018, with most of the increase occurring during 2000-2011 and the rate of increase being 17%. Finally, the proportions of the area with "poor" and "relatively poor" health remained almost unchanged. The spatial distribution of the ecosystem health in Bayanbulak is shown in Figure 7. In general, the ecosystem health levels within the property zone remained stable trend, while those within the buffer zone changed significantly. Specifically, the areas with "good" health are mainly distributed within the property zone; however, there was a decline in the health level of some sections of the northern part of the property zone in 2011. The "relatively good" health areas were mostly distributed in the buffer zone during 2000, and in the southern and eastern parts of the zone in 2011 and 2018. The areas with "ordinary" health were sporadically distributed in 2000, but their proportion increased significantly in 2011 and 2018, and they were observed in the western and northwestern parts of the study area in the buffer zone. The "poor" and "relatively poor" health areas were mainly distributed within the buffer zone while they showed dot sporadic distribution and distributed less, however, their proportion increased slightly in 2018. In summary, although the ecosystem health declined in the study area, it generally remained high. Most of the changes in the ecosystem health occurred within the buffer zone, with the health level changing from "relatively good" to "ordinary." However, the proportions of the area with "poor" and "relatively poor" health remained low. Therefore, on the whole, the ecosystem in the study area was found to be in the healthy state.

Spatial Autocorrelation Analysis
The global Moran's I index values in 2000, 2011, and 2018 were 0.762, 0.809, and 0.821, respectively, and they were all greater than 0 and showed an upward trend (Figure 8). These results were indicative of a significant positive spatial autocorrelation. Further, the spatial distributions of the ecosystem health levels showed obvious agglomeration characteristics, with the degree of agglomeration increasing over time. The LISA agglomeration chart for the local spatial autocorrelation of the ecosystem health levels is shown in Figure 9. The high-high type clusters were mainly distributed within the property zone, and their proportion increased during 2000-2018. The low-low type clusters were mainly distributed within the buffer zone, the spatial distributions during 2000 and 2011 were similar, but there was a decrease in the southern and northeastern parts in 2018, while the western and northwestern parts exhibited relatively stable distributions. In addition, the high-low and low-high type clusters were sporadically distributed.

Field Verification
The distribution of the field vegetation sampling plots in Bayanbulak is shown in Figure 1. Since none of the sampling plots were located in the areas with "poor" health, the "poor" level was not included in the field verification process. We calculated the vegetation index values for the sampling plots in each zone based on the ecosystem health levels ( Table 6). The results showed that the vegetation coverage, Simpson diversity index, Shannon-Wiener diversity index, and Margalef richness index values were consistent with the changes in the ecosystem health levels in that the higher the level of ecosystem health, the higher the values of the various vegetation indices were. The only exception was Pielou's evenness index.  Table 7 shows the results of the correlation analysis of the various vegetation indices and the ecosystem health levels performed using the Spearman rank correlation method. It was found that the r s value was greater than 0 in all the cases, indicating that the results of the field vegetation sampling process were positively correlated with the ecosystem health. The P values for the vegetation coverage, Simpson diversity index, and Shannon-Wiener diversity index were 0.019, 0.035, and 0.038, respectively. Thus, all three indices were significantly correlated with the ecosystem health. To sum up, the vegetation index values were consistent with the changes in the ecosystem health levels, and the various vegetation indices were positively correlated with the ecosystem health values. Thus, the results of RS monitoring of ecosystem health were credible.

Discussion
The health of an ecosystem is reflective of its overall state. However, absolutely healthy ecosystems do not exist, and the focus of EHA is to determine the differences in the spatial distributions in ecosystem health that occur over time [68]. The ecological effect of the land cover is a key factor in ecological protection, and has a determining effect on the health of the ecosystem [42,68,69]. The current trend in EHA research is to consider the land cover in studies on ecological health [43,44,70]. Assessing ecosystem health based on the changes in the landscape types of the ecosystem allows one to spatially understand the dynamic changes in the health levels [43]. Therefore, in this study, we focused on analysing the spatial differences and dynamic changes in ecosystem health based on the changes in the spatial distributions of the ecosystem landscape types.

Ecosystem Health Assessment model
Monitoring world heritage sites is essential for protecting them, the key to which is the construction of monitoring model [71,72]. Ecosystems are often complex, each ecosystem (such as water ecosystem and the grassland ecosystem) has different basic features [73], and the structure and function of ecological system may also be different under different space-time, therefore, many model methods have been proposed to assess the ecosystem health, and the most common methods are DMM, EMM, PSR, and VOR. Each model has different advantages and disadvantages, and they are suitable for different research areas.
DMM is based on the direct measurement and indirect calculate the value to assess ecosystem health [73], and it is often used to assess the water quality of rivers in small scale, or to assess the health of a single ecosystem [15,16]. This model could reflect the ecosystem health accurately, but it needs a lot of measured data, meaning it is unsuitable for integrated ecosystem health assessment. EMM has high requirements on the establishment of ecosystem model and the setting of model parameters [17], and it is also unsuitable for the integrated ecosystem health assessment. PSR focuses on the ecosystem state and the impact of human activities on the ecosystem and their interaction [29], it is suitable for large-medium scale region ecosystem, which is greatly affected by human activities. VOR pays more attention to the health of the ecosystem's own structure [9]. The above model often assesses health through the state of the ecosystem itself and external disturbance, and ignores the ability of the ecosystem to provide services to humans [12,27]. However, in the case of health ecosystem, it needs to meet human's reasonable requirements and keep the sustainability of itself [12]. VORS model combined the structure of the ecosystem itself and its ecosystem service to humans [12], and it is easy to describe in theory and conduct in practice [73], thus it is suitable for the study of nature reserve sites with strict restrictions on human activities.
The weights assigned to these indicators of VORS model have a determining effect on the assessment results. Since assessing ecosystem health involves humans making subjective judgments that should take into account human needs of ecosystem functional services [68], using subjective methods to determine the weights to be assigned to the indicators is one way of highlighting the relative importance of the different indicators. This is both scientifically reasonable and a widely accepted practice in EHA [12,27,42,43,56]. Therefore, we relied on expert knowledge and experience, as well as the actual conditions in Bayanbulak, which reflect its uniqueness, to determine the weights in this study. However, if this model was used in other place, it should be adjusted appropriately based on the actual conditions of study area.
The selection of data also has an important impact on the results of the model. RS data can be used to monitor and assess the spatial heterogeneity of ecosystem health at different space-time scales [50,74]. Hence, RS and GIS technologies are used widely in EHA [20,45]. RS large-scale data can quickly assess the space-time distribution and evolution of ecosystem health in heritage sites, but RS large-scale data have certain limitations in accurately reflecting the level of ecosystem health, thus, the ground-based field sampling small-scale data were used to verify the rationality of the results of RS large-scale data [65].
Bayanbulak covers the best-preserved areas of the alpine wetland ecosystem, and it has important ecological function value and provides important ecological function service value for human [40]. Thus, the ability of Bayanbulak's ecosystem to provide the value of ecological functions cannot be ignored in its ecosystem health assessment. In addition, WNH site of Bayanbulak is a nature reserve which has strict protection measures and human socioeconomic activities are restricted, and the possible negative impact of human activities is relatively small, and its direct effect can even be ignored. Therefore, we used the VORS to assess the ecosystem health of our study area from the perspective of heritage protection monitoring.
In addition, in the application of data, we used RS large-scale data and GIS analysis to assess the ecosystem health in WNH of Bayanbulak, and used the ground-based sample small-scale data to verify the rationality of the results of ecosystem health level by RS large-scale data. Therefore, our study complements the research of ecosystem health in WNH of Bayanbulak and would provide a new idea and practical experience for future monitoring and protection of WNH sites.

Ecosystem Health Assessment
The regions with "good" ecosystem health were mainly distributed within the property zone, which is a typical area with high ecosystem health, and mostly remained stable from 2000 to 2018. This was because the boundary of the study area overlaps with the NNR, and this area is strictly protected and controlled as per the laws and regulations for NNRs in China. A similar phenomenon was also found in Sanjiangyuan NNR, that areas that received strict control and protection had the highest ecosystem health level [41]. In addition, those areas contain a large number of swamp grassland and water bodies with high water content, and the land utilisation intensities were extremely low for these regions. This result is consistent with that of a study of He et al. (2019), who found that the moisture index and the land utilisation intensities were the primary parameter that determines the ecosystem health [27]. This is another reason the ecosystem health level was the highest within the property zone. Meanwhile, our field sampling data also found ecosystem health level was the highest within the property zone. In addition, we observed several rare and endangered species during the field sampling process, including Orchis umbrosa, Orchis latifolia, also suggesting that the habitat conditions in the property zone were better than those in the buffer zone.
However, the northern part of the property zone exhibited a decline in ecosystem health in 2011, and the declined areas were distributed in the restricted and exhibition zones within the property zone. A possible explanation is that those areas are close to Bayanbulak town and vulnerable to tourism and grazing activities, which would decrease the ecosystem health [75]. Jia et al. (2011) and Qian et al. (2016) also found the low ecosystem health level areas were easier to be affected by human activities in Sanjiangyuan and Dongting Lake wetland, respectively [41,76]. In addition, Bayanbulak wetland areas had begun to degenerate and shrink in 2010 [33], leading to the decline of ecosystem health in the northern part of the property zone.
The ecosystem health of the buffer zone was lower than the property zone and exhibited obvious spatial changes. We found the ecosystem health within the buffer zone in the southwestern, western, and northwestern parts were lower than the other part of the buffer zone in 2011 and 2018. This was probably because of the seasonal grazing activities in these areas from June to September. In addition, the northwestern part is adjacent to Bayinguoleng town and National Road 217. Hence, it is likely that extensive human activities in the area resulted in lower vegetation coverage [35]. In addition, we also found that the vegetation index values in the northwestern part of the study area were lower than those in the southeastern part. This was consistent with the results of RS monitoring. However, most of the buffer zone exhibited "relatively good" health and "ordinary", and there were few areas with "poor" ecosystem health, indicating that the overall health of the buffer zone was relative good. During field sampling, we found that the buffer zone contains abundance of Pedicularis kansuensis, which has low nutrient requirements and semi-parasitic characteristics [77], which would result in large-scale outbreaks and threaten the balance of the grassland ecosystem.

Limitations and Future Prospects
In this study, we established an integrated model for EHA based on the space-time changes in the ecosystem landscape types. Further, by combining RS imaging data with field sampling data, we evaluated the suitability of the proposed EHA model. However, the study has a few limitations. For instance, during the field sampling process for verifying the EHA model, the data were not collected spatially uniformly from the study area. This was owing to the inaccessibility of the wetlands. In addition, the study area shares a boundary with the Bayanbulak NNR, and the protection policy in the study area prohibited access to some of its parts. Hence, some of the sampling plots could not be accessed during the sampling process. We also found that even though the study area has strict controls and protection management measures in place, grazing and tourism activities do occur in the restricted, exhibition, and buffer zones to different degrees, which may affect the ecosystem health of Heritage. In addition, we used the RS data of 2000, 2011, and 2018 to assess the ecosystem health of Bayanbulak, and used the field data of 2018 and 2019 to verified the results of ecosystem health in 2018. As the field vegetation sampling was only conducted in 2018 and 2019, meaning we could not verify the results of ecosystem health in 2000 and 2011, which is the shortcoming of our study.
Therefore, future studies should try to ensure that the sampling plots are distributed uniformly in the study area, and the grazing and tourism activities should be accounted for within the frameworks for the assessment of heritage sites. Finally, continued efforts should be made to explore the factors driving ecological health as well as their impact mechanisms. In addition, integrated sky-space-ground monitoring systems are expected to allow for the comprehensive and focused monitoring of heritage sites at different space-time scales and provide a scientific basis for the establishment of suitable protection measures.

Conclusions
In this study, we established an indicator-based framework for the EHA of WNH sites and evaluated it using Bayanbulak in China as a test case for the periods of 2000, 2011, and 2018. We analysed the space-time evolution characteristics and spatial agglomeration effects of the ecosystem health levels based on the changes in the ecosystem landscape types from 2000 to 2018. The primary results of this study can be summarised as follows:

•
While changes were observed in the area proportions of the various ecosystem landscape types from 2000 to 2018, the overall landscape structure did not change. The fragmentation of the landscape in the study area increased in general during the three periods. However, the degree of fragmentation remained small. • All the EHA indicators showed a decline, with the decline in ES being the most significant. Thus, the overall ecosystem health in the study area also declined. However, the ecosystem health within the property zone remained high and essentially unchanged. Further, the area proportions of poor health were extremely low and were mostly distributed within the buffer zone. Therefore, in general, the ecosystem of the study area was in a healthy state.

•
The spatial distribution of the ecosystem health exhibited obvious agglomeration characteristics, with the degree of agglomeration increasing over time. With respect to the local spatial autocorrelation, the high-high clusters were mainly distributed within the property zone, and the degree of agglomeration enhanced over time. On the other hand, the low-low clusters were mainly distributed in the buffer zone, and the degree weaken over time.

•
The results of EHA obtained through RS monitoring were positively correlated with those of the field sampling of the vegetation in the study area. The areas with high levels of ecosystem health also exhibited high sampling vegetation index values. This confirmed the suitability of using RS imaging for monitoring ecosystem health.
We diagnosed the ecosystem health level in WNH of Bayanbulak and analyzed the characteristics of its space-time evolution, which helps local heritage managers to formulate precise protection measures for different regions with different ecosystem health levels, and provided guiding significance for the protection of Bayanbulak. In addition, we established an EHA framework in the WNH sites, which could be applied to other nature reserve sites for their better protection in the future. However, the RS data used to assess the ecosystem health were large-scale data, which may affect the accuracy of results, so we used small-scale field data to verify its rationality. However, due to the lack of field sampling data in 2000 and 2011, the accuracy of EHA was not fully verify. In the future, we could increase the accuracy of EHA framework in WNH by using more ground-based sample small-scale data and unmanned Aerial Vehicle medium-scale data.