Ecological Environment Assessment in World Natural Heritage Site Based on Remote-Sensing Data. A Case Study from the Bayinbuluke

Ecological environment assessment would be helpful for a rapid and systematic understanding of ecological status and would contribute to formulate appropriate strategies for the sustainability of heritage sites. A procedure based on spatial principle component analysis was employed to measure the ecological status in Bayinbuluke; exploratory spatial data analysis and geo-detector model were introduced to assess the spatio-temporal distribution characteristics and detect the driving factors of the ecological environment. Five results are presented: (1) During 2007–2018, the average values of moisture, greenness, and heat increased by 51.72%, 23.10%, and 4.99% respectively, and the average values of dryness decreased by 56.70%. However, the fluctuation of each indicator increased. (2) The ecological environment of Bayinbuluke was improved from 2007 to 2018, and presented a distribution pattern that the heritage site was better than the buffer zone, and the southeast area was better than the northwest area. (3) The ecological environment presented a significant spatial clustering characteristic, and four types of spatial associations were proposed for assessing spatial dependence among the samples. (4) Elevation, protection partition, temperature, river, road, tourism, precipitation, community resident, and slope were statistically significant with respect to the changes in ecological status, and the interaction of any two factors was higher than the effect of one factor alone. (5) The remote-sensing ecological index (RSEI) could reflect the vegetation growth to a certain extent, but has limited ability to respond to species structure. Overall, the framework presented in this paper realized a visual and measurable approach for a detailed monitoring of the ecological environment and provided valuable information for the protection and management of heritage sites.


Introduction
Natural World Heritage Sites (NWHSs) represent areas of outstanding universal value from the view of aesthetic importance, geology, ecosystems or biodiversity [1][2][3]. A global consensus has been reached that NWHSs are significant components in global ecosystem protection which is so exceptional as to transcend national boundaries and to be of common importance for present and future generations of all humanity [1]. However, in recent years, NWHSs have been experiencing unprecedented pressures due to climate change, natural disaster, and human activities [4][5][6][7][8][9].

Study Area
The Bayinbuluke is located in Hejing County, Xinjiang, China (central coordinates: N42 • 47 53 , E84 • 09 50 ), with a total area of 109,448 ha and a buffer zone of 80,090 ha ( Figure 1). The study area is situated in the intermontane Youerdusi Basin, and the Kaidu River originating from Aierbin Mountain crosses the center of the basin, forming a broad floodplain and swamp grasslands. In the arid continental climate zone, the study area has a short summer, a long and freezing winter, and a very short frost-free period all year round. The annual average temperature is −4.6 • C, and the annual average precipitation is 276 mm, mainly concentrated from June to August [32]. Because of frosty weather and frozen soil, there are hardly any large trees in this region, the main vegetation types are Alpine steppe and Alpine marsh meadow. Bayinbuluke is the best habitat for breeding and a summer habitat for wild birds. There are abundant species of birds in Bayinbuluke with in total 119 species of all kinds, such as cranes, egrets, golden eagles, bar-headed geese, grey geese, whooper swans, tundra swans and mute swans. The Bayinbuluke is located in Hejing County, Xinjiang, China (central coordinates: N42°47'53", E84°09'50"), with a total area of 109,448 ha and a buffer zone of 80,090 ha ( Figure. 1). The study area is situated in the intermontane Youerdusi Basin, and the Kaidu River originating from Aierbin Mountain crosses the center of the basin, forming a broad floodplain and swamp grasslands. In the arid continental climate zone, the study area has a short summer, a long and freezing winter, and a very short frost-free period all year round. The annual average temperature is −4.6 ℃, and the annual average precipitation is 276 mm, mainly concentrated from June to August [32]. Because of frosty weather and frozen soil, there are hardly any large trees in this region, the main vegetation types are Alpine steppe and Alpine marsh meadow. Bayinbuluke is the best habitat for breeding and a summer habitat for wild birds. There are abundant species of birds in Bayinbuluke with in total 119 species of all kinds, such as cranes, egrets, golden eagles, bar-headed geese, grey geese, whooper swans, tundra swans and mute swans.
The heritage values of Bayinbuluke can be summarized as follows: (1) Bayinbuluke is the best representative of the large inter-montane basins of the Tianshan Mountains. (2) Bayinbuluke is the typical representative of an alpine wetland ecosystem in the arid temperate zone. (3) Bayinbuluke is the best representative of a beautiful landscape of bending rivers and marshes of the Tianshan Mountains. (4) Bayinbuluke is China's largest breeding ground for swans, as well as the northern hemisphere's most southern limit for the breeding of wild swans.

Data Source and Pre-Processing
To evaluate the ecological environment quality in Bayinbuluke, remote-sensing data, spatial data, the material declarations for Xinjiang Tianshan NWHS, and filed survey data were collected. The heritage values of Bayinbuluke can be summarized as follows: (1) Bayinbuluke is the best representative of the large inter-montane basins of the Tianshan Mountains. (2) Bayinbuluke is the typical representative of an alpine wetland ecosystem in the arid temperate zone. (3) Bayinbuluke is the best representative of a beautiful landscape of bending rivers and marshes of the Tianshan Mountains. (4) Bayinbuluke is China's largest breeding ground for swans, as well as the northern hemisphere's most southern limit for the breeding of wild swans.

Data Source and Pre-Processing
To evaluate the ecological environment quality in Bayinbuluke, remote-sensing data, spatial data, the material declarations for Xinjiang Tianshan NWHS, and filed survey data were collected. Each thematic layer was processed with the same resolution (30 meter), projected coordinate system (WGS_1984_UTM_45N) and spatial boundary. (1) Remote-sensing data. Landsat images (column number:145; row number:30) in 2007 and 2018 were downloaded from the United States Geological Survey (https://www.usgs.gov/). The images were preprocessed with the ENVI 5.3 image processing software. Radiometric correction converted the digital number of the raw images to reflectance at sensor, and then, FLAASH atmospheric correction was employed to reduce the deviations caused by light and atmosphere. (2) Spatial data. Digital elevation model (DEM) data and thematic maps about temperature and precipitation were obtained from the Resources and Environment Data Cloud Platform (http://www.resdc.cn/). Altitude and slope were retrieved from DEM data.
(3) The material declarations for Xinjiang Tianshan NWHS. The material contained the Bayinbuluke nomination text and spatial information on the river, road, tourism attractions, and community settlements of the study area. (4) To verify whether the results were in accordance with reality, a field survey was conducted from July 18, 2018 to July 23, 2018. There were 31 samples taken close to the villages and tourism attractions were selected due to budget constraints and poor road accessibility. Grass vegetation samples were collected in the main transects of 1 m × 1 m, and then, plant species, quantity, height, coverage, and aboveground biomass were measured. The Margalef richness index (Rm), Shannon-Wiener diversity index (H'), Simpson diversity index (D), Pielou evenness index (JSW), and Gini uniformity index (Jgi) of each sample was calculated with regard to species diversity index [33,34].

Ecological Indicators Extracted
Land surface information extracted from remote sensing images are widely used to evaluate ecological conditions because it is an effective describing way to ecological quality. The remote-sensing ecological index (RSEI) is constructed to assess ecological environment [23,35].
• Land surface moisture Tasseled cap transformation has been extensively applied in ecological assessment studies [36], due to the components of Tasseled cap having a direct impact on the physical parameters of the earth's surface [37]. As a wetness component, the land surface moisture (LSM) was calculated based on the reflectance of bands. The calculation of LSM is expressed as follows: is the planetary reflectance of blue, green, red, near-infrared and mid-infrared band of the Landsat 5 TM image and Landsat 8 OLI image respectively.

• Greenness
Vegetation growth describes the confrontation between vegetation and natural environment [38]. Normalized Difference Vegetation Index (NDVI) was employed to measure vegetation coverage and growth with near-infrared band and red band for the capacity to represent greenness. [39,40]. NDVI can be calculated as follows: where b Red and b NIR represent the planetary reflectance of red band and near-infrared band respectively. • Dryness Normalized Difference Imperviousness and Soil Index (NDISI) can effectively distinguish impervious and soil features without covering water before processing impervious surface information [41]. The Impervious Built-up Index (IBI) has been commonly used to map built-up lands accurately. In addition to the built-up lands, patches of bare land or sparsely vegetated ground occurred in the deforested or abandoned locations across the study area. For this reason, Soil Index (SI) was also employed to represent these bare areas. NDISI was composed of IBI and SI, as proposed here: is the planetary reflectance of blue, green, red, near-infrared and mid-infrared band of the Landsat 5 TM image and Landsat 8 OLI image respectively.
• Land surface temperatures A unique feature of thermal infrared band in remote sensing data is the capability to measure land surface temperatures (LST) [42]. All satellite-based studies of heat island have verified a close relationship between the satellite-derived surface temperature and land surface temperature [43,44]. The calculations of LST were expressed as follows [45]: where L λ is the at-satellite spectral radiance values of the thermal bands, i.e., Band 6 of Landsat TM and Band 10 of Landsat OLI; gain and bias is the gain value and offset value of the corresponding band, both of the values can be queried from the header file of corresponding image. T b is the at-satellite brightness temperature, K 1 and K 2 are the band-specific thermal conversion constants of thermal bands, at Band 6 of Landsat TM, K 1 = 607.76 W/ m 2 · sr· µm , K 2 = 1260.56K; at Band 10 of Landsat OLI, K 1 = 774.89 W/ m 2 · sr· µm , K 2 = 1321.08 K. λ is the wavelength of the thermal infrared band; ρ = 1.4380 * 10 4 µm; ε is the surface specific emissivity.

Construction of Remote-Sensing Ecological Index (RSEI)
The initial data was standardized in a range from 0~1 to remove the influence of different measurement units. For the index which has positive correlation with RSEI, the equation is expressed as follows: While, for the index which has negative correlation, the equation is: where, X i is the standardized value of variable i, x i is the measured value of variable i. Principal component analysis is skilled in reducing the dimensionality of a dataset by converting a set of observed correlated variables into a set of linearly uncorrelated variables through orthogonal transformation [46]. Moreover, spatial principal component analysis does not rely on the prior knowledge or experience of researchers, which can reduce subjective influence to some extent [47]. Remote-sensing data and geographic information system provide an available technique in integrating and analyzing a series of spatial data. Regarding the difficulty in quantifying spatial and temporal variations, the spatial principal component analysis, integrating the geographic information system and principal component analysis [48], was employed here. The principal components are selected based on the fact that the first principal component represents the greatest amount of variance in the data. If the accumulated variance represents over 85% of the total, the remaining components can be ignored. The spatial principal component analysis provided objective results in the process of evaluating the ecological environment. A synthetic index (i.e., RSEI) that will allow for a quick and quantitative assessment of ecological environment is calculated as follows: where i is the quantity of principal component (PC) that remained and r i is the contribution ratio of PC i . The contribution ratio r i is calculated as follows: where p i represents the contribution ratio of principal component p i , and n is the significant number of principal components that remain. Based on the result of spatial principal component analysis (Table 1), the formulas of RSEI in 2007 and 2018 are expressed as: The distribution maps of ecological status in 2007 and 2018 were drawn by the raster calculation tool in Arcgis10.5 software. Referring to the characteristics of the Bayinbuluke, the values were divided into five grades which represent very bad, bad, moderate, good, excellent, respectively, i.e., RSEI < 0.8: Grade I (very bad); 0.8 ≤ RSEI < 0.9: Grade II (bad); 0.9 ≤ RSEI < 1.0: Grade III (moderate); 1.0 ≤ RSEI < 1.1: Grade IV (good); RSEI ≥1.1: Grade V (excellent).
The Remote-Sensing Ecological Body Index (RSEBI) was employed to estimate the overall differences in ecological status, and it is calculated as follows: where P i is the ecological environment grade; 1 to 5 are assigned to Grade I to Grade V respectively; A i denotes the area of grade i; S is the total area of Bayinbuluke.

Exploratory Spatial Data Analysis
Exploratory spatial data analysis is a collection of methods used to describe and visualize the spatial distribution of an attribute. It is also used to suggest spatial regimes or other forms of spatial heterogeneity [49] by identifying atypical locations or spatial outliers and discovering the patterns of spatial association. In practical application, global spatial autocorrelation and local indicator of spatial association (LISA) are often adopted to investigate the spatial characteristics of observations.

• Global Spatial Autocorrelation
Global spatial autocorrelation is employed to measure if and how much the attribute is autocorrelated in the whole region. Global Moran's I statistics is defined by the following [50,51]: x i ; n is the number of spatial units, x i and x j are the observations of spatial units i and j, respectively; w ij is an element of the spatial weight matrix W which describes the spatial arrangement of all the spatial units in the sample, if spatial units i and j were adjacent, w ij is 1, otherwise, w ij is 0. The values of Global Moran's I range from −1 to 1.
• Local indicator of spatial association (LISA) LISA is the most commonly used indicator among the local indicators for spatial association [52]. It can be used to determine the correlation degree of an attribute between a unit and its adjacent units throughout the region. LISA is useful to identify local spatial cluster patterns and spatial outliers [53], and it was defined as:

Geo-Detector Model
The geo-detector model is new statistical method to detect spatial stratified differentiation among the geographical elements and reveal the driving factors behind it. The study area is characterized by spatial stratified heterogeneity if the sum of the variance of subareas is less than the regional total variance; and if the spatial distribution of the two variables tends to be consistent, there is statistical correlation between them [54]. The Q-statistic in the geo-detector has already been applied in many fields of natural and social sciences which can be used to measure spatial stratified heterogeneity, detect explanatory factors, and analyze the interactive relationship between variables. In this paper, factor detection and interaction detection were employed to analyze the driving force and the interaction between multiple elements. The level of which geographical element X explained the spatial variation of attribute Y was detected by the following formula: where q is the explaining power of a variable on a spatial attribute; t is the categories or partitions of variables; N h is the quantity of sample units in subfields; n is the quantity of sample units in the whole area; L is the quantity of subfields; σ 2 is the variance of a single variable of the entire region; and σ 2 h is the variance of the subfield. Interaction judgement has the ability to detect whether the various variables formed a combined effect on the ecological environment, or whether the influence between different variables was independent. The basis for interaction type is shown in Table 2. Table 2. Types of interaction between two covariates.

Variables Standard of Classification
Elevation (×1) Extracted from digital elevation model (DEM) data and classed into five categories with natural breaks. Slope (×2) Extracted from DEM data and classed into five categories with natural breaks. Precipitation (×3) Classed into five categories with natural breaks. Temperature (×4) Classed into five categories with natural breaks.

River (×5)
Kernel density of river is created, and classed into five categories with natural breaks.

Road (×6)
Kernel density of road is created, and classed into five categories with natural breaks.
Community resident (×7) Kernel density of town and grazing point is created, and classed into five categories with natural breaks.
Tourism (×8) Kernel density of tourism infrastructures is created, and classed into five categories with natural breaks.
Protection partition (×9) In 2007, Bayinbuluke was divided into experimental area, core area and buffer zone according to national nature reserve requirements. In 2018, Bayinbuluke was divided into construction prohibited district, construction limited district, exhibition district and buffer zone according to the management plan of Xinjiang Tianshan.

Spatial Econometric Model
After determining the spatial correlation of the regional unit attributes, the spatial interaction between different areas should be introduced into the model as a variable to construct a spatial econometric model, which mainly includes three types of model: ordinary least squares (OLS), spatial lag model (SLM) and a spatial error model (SEM). In the spatial correlation test, if the spatial lag model Lagrange multiplier test statistic LMLAG and spatial error LMERROR are not significant, then the OLS regression were selected. If LMLAG is more significant than the spatial error LMERROR, the SLM is selected, otherwise, the SEM is selected [47,55].
The SLM can be calculated as follows: where, y is the dependent variable matrix; X is the argument matrix; the parameter β reflects the influence of the independent variable on the dependent variable; W is a spatial weights matrix; and ε represents the residual items.
The SEM can be calculated as follows: where, λ is the coefficient of autocorrelation for the SEM.

Characteristics of Land Surface Moisture (LSM), Normalized Difference Vegetation Index (NDVI), Normalized Difference Imperviousness and Soil Index (NDISI) and Land Surface Temperatures (LST)
From the spatial distribution of ecological factors (see in Figure 2), LSM and NDVI were higher in the eastern and lower in the western areas, which were closely related to the distribution of the river system. By contrast, the values of NDISI and LST were higher in the western part and lower in the eastern part. In 2007-2018, the moisture and temperature of land surface were decreased in the low vegetation coverage area in west and north, and the improvement of vegetation cover in southeast and the middle swamp and wetland area was obvious. The statistics of four indicators are shown in Table 4. The average value of LSM, NDVI, LST increased during the period of 2007-2018, while NDISI decreased. Among them, the increase in LSM was large, with a rate of 51.72%, followed by LST, with a rate of 23.10%, and NDVI, with a rate of 0.81%. The NDISI was on the decrease, and its mean value decreased from −0.0866 in 2007 to −0.1357 in 2018, indicating a significant decrease in built-up land and bare land in the study area. The standard deviation of LSM, NDVI, NDISI and LST increased, indicating that the extreme differences in moisture, greenness, dryness and heat tend to be obvious.

Characteristics of Ecological Environment
Redundancy of data occurs frequently and should be considered when processing data [56]. The variance inflation factor (VIF) and tolerance index (TOL) were employed to measure the potentially redundant features for the purpose of effectiveness of dimensionality reduction. If VIF > 10 or TOL < 0.1, this indicated that the disturbance of multicollinearity in the dataset was serious [57]. VIF and TOL were computed with the sample points related to RSEI and principal components in 2007 and 2018 (see in Table 5). The result showed that VIF < 10 and TOL > 0.1, both of which illustrated that the redundancy among the principal components were acceptable, and the selection of indicators and methods were desirable. Table 5. Results of variance inflation factor (VIF) and tolerance index (TOL). The statistics of four indicators are shown in Table 4. The average value of LSM, NDVI, LST increased during the period of 2007-2018, while NDISI decreased. Among them, the increase in LSM was large, with a rate of 51.72%, followed by LST, with a rate of 23.10%, and NDVI, with a rate of 0.81%. The NDISI was on the decrease, and its mean value decreased from −0.0866 in 2007 to −0.1357 in 2018, indicating a significant decrease in built-up land and bare land in the study area. The standard deviation of LSM, NDVI, NDISI and LST increased, indicating that the extreme differences in moisture, greenness, dryness and heat tend to be obvious.

Characteristics of Ecological Environment
Redundancy of data occurs frequently and should be considered when processing data [56]. The variance inflation factor (VIF) and tolerance index (TOL) were employed to measure the potentially redundant features for the purpose of effectiveness of dimensionality reduction. If VIF > 10 or TOL < 0.1, this indicated that the disturbance of multicollinearity in the dataset was serious [57]. VIF and TOL were computed with the sample points related to RSEI and principal components in 2007 and 2018 (see in Table 5). The result showed that VIF < 10 and TOL > 0.1, both of which illustrated that the redundancy among the principal components were acceptable, and the selection of indicators and methods were desirable. According the principal component analysis method, the values of RSEI were calculated and divided into five grades. As seen in Figure 3, ecological environment presented a distribution pattern that the NWHS was better than the buffer zone, and the southeast area was better than the northwest area. In 2007, the ecological environment grade was dominated by Grade III, mainly distributed in the heritage site and the northeast of the buffer zone. The areas with Grade IV and Grade V were embedded in the Grade III, mainly distributed along the rivers. The areas with Grade II were mainly distributed in the buffer zone, and Grade I areas were mainly concentrated in the northern of the buffer zone, where the density of rivers was scarce, and the main landscape was low coverage grassland. In 2018, the ecological environment grade was dominated by Grade IV, mainly distributed in the heritage site and the southeast of the buffer zone. Most of the areas with Grade III were replaced by the areas with Grade IV and Grade V. The areas with Grade I and Grade II were mainly concentrated in the west of the buffer zone, and their area were significantly reduced. The ecological environment of Bayinbuluke has been improved as a whole, the RSEBI increased from 2.7017 in 2007 to 3.5225 in 2018. According the principal component analysis method, the values of RSEI were calculated and divided into five grades. As seen in Figure 3, ecological environment presented a distribution pattern that the NWHS was better than the buffer zone, and the southeast area was better than the northwest area. In 2007, the ecological environment grade was dominated by Grade III, mainly distributed in the heritage site and the northeast of the buffer zone. The areas with Grade IV and Grade V were embedded in the Grade III, mainly distributed along the rivers. The areas with Grade II were mainly distributed in the buffer zone, and Grade I areas were mainly concentrated in the northern of the buffer zone, where the density of rivers was scarce, and the main landscape was low coverage grassland. In 2018, the ecological environment grade was dominated by Grade IV, mainly distributed in the heritage site and the southeast of the buffer zone. Most of the areas with Grade III were replaced by the areas with Grade IV and Grade V. The areas with Grade I and Grade II were mainly concentrated in the west of the buffer zone, and their area were significantly reduced. The ecological environment of Bayinbuluke has been improved as a whole, the RSEBI increased from 2.7017 in 2007 to 3.5225 in 2018. The area of each ecological status grade was calculated and shown in Figure 4. In 2007, the areas with Grade I to Grade V were 217.30km 2 , 568.21km 2 , 698.22km 2 , 385.84km 2 , and 25.80km 2 , respectively. From 2007 to 2018, the areas with Grade IV and Grade V increased to 673.11km 2 and 444.09km 2 , respectively, while the areas with Grade I, Grade II, and Grade III decreased to 128.85km 2 , 316.27km 2 , and 330.06km 2 , respectively. Area conversion mainly occurred between adjacent levels, and presented a conversion from the low level to the high level. The main types of shift included Grade V to Grade IV, Grade IV to Grade III, Grade III to Grade II, Grade III to Grade I, Grade II to Grade I. The area of each ecological status grade was calculated and shown in Figure 4. In 2007, the areas with Grade I to Grade V were 217.30 km 2 , 568.21 km 2 , 698.22 km 2 , 385.84 km 2 , and 25.80 km 2 , respectively. From 2007 to 2018, the areas with Grade IV and Grade V increased to 673.11 km 2 and 444.09 km 2 , respectively, while the areas with Grade I, Grade II, and Grade III decreased to 128.85 km 2 , 316.27 km 2 , and 330.06 km 2 , respectively. Area conversion mainly occurred between adjacent levels, and presented a conversion from the low level to the high level. The main types of shift included Grade V to Grade IV, Grade IV to Grade III, Grade III to Grade II, Grade III to Grade I, Grade II to Grade I.

Spatial Clustering of Ecological Environment
The RSEI map in 2007-2018 were used to examine the Global Moran's I, which can describe the overall correlation. According to Figure 5, the Global Moran's I passed the significance test (at 0.02 significance level) in all years, implying that the ecological environment quality in Bayinbuluke had significant spatial autocorrelation characteristics. From the perspective of the evolution trend, the Global Moran's I from 2007 to 2018 showed the characteristics of increasing, the spatial clustering degree of ecological environment in the research area was continuously strengthened, and the clustering level was enhanced. The area of each ecological status grade was calculated and shown in Figure 4. In 2007, the areas with Grade I to Grade V were 217.30km 2 , 568.21km 2 , 698.22km 2 , 385.84km 2 , and 25.80km 2 , respectively. From 2007 to 2018, the areas with Grade IV and Grade V increased to 673.11km 2 and 444.09km 2 , respectively, while the areas with Grade I, Grade II, and Grade III decreased to 128.85km 2 , 316.27km 2 , and 330.06km 2 , respectively. Area conversion mainly occurred between adjacent levels, and presented a conversion from the low level to the high level . The main types of shift included  Grade V to Grade IV, Grade IV to Grade III, Grade III to Grade II, Grade III to Grade I, Grade II to  Grade I   The RSEI map in 2007-2018 were used to examine the Global Moran's I, which can describe the overall correlation. According to Figure 5, the Global Moran's I passed the significance test (at 0.02 significance level) in all years, implying that the ecological environment quality in Bayinbuluke had significant spatial autocorrelation characteristics. From the perspective of the evolution trend, the Global Moran's I from 2007 to 2018 showed the characteristics of increasing, the spatial clustering degree of ecological environment in the research area was continuously strengthened, and the clustering level was enhanced. The Local Moran's I index was employed to assess the spatial dependence among the samples. Four types of spatial associations were proposed as follows: High-High Cluster type, High-Low Outlier type, Low-High Outlier type, and Low-Low Cluster type. High-High and Low-Low Cluster correspond to positive spatial autocorrelation, while High-Low and Low-High Outlier correspond to negative spatial autocorrelation. The characteristics of local spatial agglomeration were summarized: ① High-High Cluster type. The spatial difference of the High-High Clustering type was small. The values of samples and their adjacent samples were high, and they showed a significant positive correlation. Most of the High-High Clustering samples were distributed in the middle of the study area in the areas of Aerle, Bayintala, Bawuerken, Haerwusu, Wulanwusu, Yikeyuwulezen, and showed a blocky distribution. ② High-Low Outlier type. The sample itself has higher ecological quality, while its adjacent samples have lower ecological quality, showing a negative correlation of "high itself, low surrounding". High-Low Outlier samples were scattered in the fringes of the study area, and showed a pointlike distribution structure. ③ Low-High Outlier type. The value of the sample was low, and the values of the adjacent samples were high. Low-High Outlier type showed a negative correlation of "low itself, high surrounding", and they mainly distributed along the river. ④ Low-Low Cluster type. The spatial difference of Low-Low Clustering type was small. The samples and their adjacent samples were both in low ecological quality, and showed a significant positive correlation. The number of "Low-Low" Cluster samples gradually decreased, and they mainly distributed in Aolunbuluke, Wulalianyinge, Halahatecahan, Haersala, Aerxiate, and showed a blocky distribution. The Local Moran's I index was employed to assess the spatial dependence among the samples. Four types of spatial associations were proposed as follows (see Figure 6): High-High Cluster type, High-Low Outlier type, Low-High Outlier type, and Low-Low Cluster type. High-High and Low-Low Cluster correspond to positive spatial autocorrelation, while High-Low and Low-High Outlier correspond to negative spatial autocorrelation. The characteristics of local spatial agglomeration were summarized: 1 High-High Cluster type. The spatial difference of the High-High Clustering type was small. The values of samples and their adjacent samples were high, and they showed a significant positive correlation. Most of the High-High Clustering samples were distributed in the middle of the study area in the areas of Aerle, Bayintala, Bawuerken, Haerwusu, Wulanwusu, Yikeyuwulezen, and showed a blocky distribution. 2 High-Low Outlier type. The sample itself has higher ecological quality, while its adjacent samples have lower ecological quality, showing a negative correlation of "high itself, low surrounding". High-Low Outlier samples were scattered in the fringes of the study area, and showed a point-like distribution structure. 4 Low-Low Cluster type. The spatial difference of Low-Low Clustering type was small. The samples and their adjacent samples were both in low ecological quality, and showed a significant positive correlation. The number of "Low-Low" Cluster samples gradually decreased, and they mainly distributed in Aolunbuluke, Wulalianyinge, Halahatecahan, Haersala, Aerxiate, and showed a blocky distribution.
All factors were statistically significant with respect to the changes in the ecological status. Precipitation, temperature, road, community resident, tourism and protection partition were selected to determine the interaction effect on the changes in ecological environment for the larger variations of their absolute change and relative change. The interaction detection result was higher than the effect of one variable alone, and the interaction of any two factors showed non-linear strengthening or mutual strengthening ( Table 7). The greatest interaction was the non-linear enhancement of precipitation ∩ protection partition, which meant that the difference of ecological status in different precipitation under the same protection partition or different protection partition at the same precipitation was the largest. In 2007, the interaction of precipitation ∩ protection partition (0.4235), temperature ∩ protection partition (0.3829), and community resident ∩ protection partition (0.3625) were the significant control factors for the ecological condition in the Bayinbuluke. In 2018, the top three significant control factors were the interaction of precipitation ∩ protection partition (0.4812),

Relationship between Detected Factors and Ecological Environment
According to the analysis of the geographical detector (Table 6)  All factors were statistically significant with respect to the changes in the ecological status. Precipitation, temperature, road, community resident, tourism and protection partition were selected to determine the interaction effect on the changes in ecological environment for the larger variations of their absolute change and relative change. The interaction detection result was higher than the effect of one variable alone, and the interaction of any two factors showed non-linear strengthening or mutual strengthening ( Table 7). The greatest interaction was the non-linear enhancement of precipitation ∩ protection partition, which meant that the difference of ecological status in different precipitation under the same protection partition or different protection partition at the same precipitation was the largest. In 2007, the interaction of precipitation ∩ protection partition (0.4235), temperature ∩ protection partition (0.3829), and community resident ∩ protection partition (0.3625) were the significant control factors for the ecological condition in the Bayinbuluke. In 2018, the top three significant control factors were the interaction of precipitation ∩ protection partition (0.4812), tourism ∩ protection partition (0.3531), and temperature ∩ protection partition (0.3380).

The Correlation between RSEI and Survey Data
In case there is spatial heterogeneity between samples, the spatial regression model should be used to avoid the estimation error. Therefore, we used the spatial regression model to discuss the spatial heteroskedasticity and spatial dependence of survey data for the analysis of remote-sensing results. Then, an OLS model, a SLM, and a SEM were established. To choose the best model fitting the data, the regression models were evaluated by the results from the Lagrange Multiplier tests [58]. In this test, both LMlag (LMlag = 4.3842, P = 0.0363) and LMerr ((LMerr = 4.6380, P = 0.0313) were significant for RSEI. Then, the robust model forms (i.e., RLMerr and RLMlag) was considered. The statistic of RLMerr (0.8716) was higher than RLMlag (0.6178), which indicates the SEM was a better alternative. Moreover, the Akaike information criterion (AIC) of SEM (−88.9148) was lower than OLS models (−80.1745) and SLM (−81.8813). This criterion is based on maximum likelihood function, and lower values denote best model fit [58,59]. Table 8 presented the result of standardized regression coefficients from the SEM. In this model, aboveground biomass, coverage, and Shannon diversity index presented significant positive association with the remote-sensing result, while the Simpson diversity index had a significant and negative coefficient. In addition, Margalef richness index, Pielou evenness index, and Gini evenness index were not significant among the independent variables. The estimated parameters and their significance in the regression model revealed that the remote-sensing result could reflect the basic situation of vegetation growth to a certain extent, but had limited ability to respond to species structure such as evenness or richness.

Driving Mechanism of Ecological Environment
Altitude and precipitation are the mainly natural factors affecting the ecological environment of Bayinbuluke. Elevation itself has a correlation with climate and topographic factors, and affects the distribution of vegetation, river and landscape patterns on a large scale [60,61]. In the material declarations for Xinjiang Tianshan NWHS, the natural vegetation in Bayinbuluke can be divided into Alpine steppe and Alpine marsh meadow from the high to the low altitude. Alpine steppe includes Stipa purpurea steppe, Stipa grandis steppe and miscellaneous grass prairie. This kind of vegetation grows at an elevation of 2500-2600 m, and the community is composed of cold-tolerant plants such as Stipa subsessiliflora, S.purpurea, S.krylovii, Festuca pseudovina, Gentiana decumbens, potentilla bifurca, Agropyron cristatum, poa sp. and Schultzia sp. On the slopes and pediment of the mountains, there is widespread alpine Stipa purpurea steppe; on the wet basin floor, there is marsh meadow and marsh vegetation. Alpine marsh meadow mainly includes Carex meadow which is located in the basin bottom. Herb species growing in marsh meadow account for more than 60% of all plant species in the heritage site, such as Carex uesicaria and C.microglochin, as well as wet and mesic grass such as Potamogeton lucens, Utricularia sp., Hippuris vulgaris, Trigloch in palustre, Alopecurus arundinaceus and Deschampsia caespitosa. There is also an abundant distribution of submerged plants and emergent plants in the lakes and rivers, as well as Utricularia, Halerpestes cymbalaris and Potamogeton malaianus in the water. Precipitation of the growing season (May-September) in 2018 increased by 9.88% compared with 2007, and the RSEBI of 2018 is greater than 2007. The impact of precipitation on the ecological environment in our study supports the findings of some previous studies. For instance, Zhao et al. [62] studied the relationship between interannual vegetation change and climate factors, and found that increasing precipitation induced an increase in NDVI in general. Mo et al. [63] concluded that most of the area in the arid basin shows a strong positive correlation between precipitation and NDVI at different temporal and spatial scales.
Anthropogenic factors include conservation management and tourism activities. The management plan of Xinjiang Tianshan divided Bayinbuluke into construction prohibited district, construction limited district, exhibition district and buffer zone. Then, corresponding protection requirements for each district was proposed. Construction prohibited district and construction limited district are the main distribution district of high RSEI, and human activities such as tourism and grazing are prohibited here. As the hotpots of biodiversity and major habitat for species, these districts are typical areas that reflect the biological and ecology value of Bayinbuluke. With the improvement of protection management requirements, tourism activities have gradually concentrated in the exhibition district, and the construction of tourism infrastructures and roads have had a certain negative impact on the ecological environment [64][65][66]. In addition, grazing activities are allowed in buffer zone. The adaptability of vegetation species to grazing disturbances is limited, the carrying capacity and the rehabilitation of grassland depend on grazing intensity, grazing strategies, and different livestock [67][68][69]. Therefore, grazing also has a complex impact on the ecological quality of the buffer zone.

Disadvantage and Further Advancement
Indicators extracted from remote-sensing image allowed us to achieve a rapid and objective evaluation and a more detailed understanding of the ecology environment of heritage sites for the superiority of spatial visualization. Considering the important values of the study area in terms of biodiversity and ecosystem, the field survey data verified RSEI from the aspects of meadow coverage, aboveground biomass, and plant community diversity. The evaluation results obtained by remote-sensing data have a positive relationship with vegetation growth, however, the relationship of vegetation community structure and remote-sensing result was not obvious, which may be detrimental to the diagnosis of degradation caused by changes in vegetation types, such as species invasion and excess population. More samples should be surveyed to explore the relationship between different levels of ecological environment and plant community diversity. With the development of monitoring techniques and tools, the equipment such as remote-sensing satellites, unmanned aerial vehicle, and infrared camera provides multi-scale and multi-dimensional monitoring data for heritage site monitoring. A remote-sensing-unmanned aerial vehicle-ground survey synergy monitoring system should be taken into account in future research. In addition, the outstanding universal value of the site and the management plan of heritage sites impose requirements on the function of different regions. For instance, the function of a construction-prohibited district is the main area for the protection of outstanding universal value, while exhibition district provides services to support tourism, local residents, and management. Therefore, protection management requirements and main functions of different regions should be taken into account for a more pertinent evaluation.

Conclusions
Ecological status assessment is significant to improve the management of heritage sites and ensure their sustainability. This study attempted to assess the ecological environment with remote-sensing data in a raster scale, which would be helpful for a rapid and systematic understanding of NWHSs. The results showed that the moisture, greenness, and heat increased, and the dryness decreased during 2007-2018 in Bayinbuluke, while the extreme differences of the four indicators tended to be apparent. The ecological status of Bayinbuluke had been improved as a whole, and the ecological environment quality presented a distribution pattern that the NWHS was better than the buffer zone, the southeast area was better than the northwest area, and the ecological grades presented a conservation from the low level to the high level. There was a significant spatial autocorrelation characteristics of ecological environment quality in Bayinbuluke, and the clustering degree was continuously strengthened. Elevation, protection partition, temperature, river, road, tourism, precipitation, community resident, and slope were statistically significant with respect to the ecological status. The interaction of any two factors was higher than the effect of one variable alone, and showed non-linear strengthening or mutual strengthening. The remote-sensing result reflected the basic situation of vegetation growth to a certain extent, but had limited ability to response species structure such as evenness or richness. The framework presented in this paper provides a visual and measurable approach for a detailed understanding of ecological environment, and the method can be employed to other heritage sites with suitable datasets, particularly in the sites representing significant on-going ecological/biological processes and development of ecosystems. Furthermore, multi-dimensional data and protection management requirements should be taken into account for a more pertinent evaluation.