Analysis of Eco-Environmental Quality of an Urban Forest Park Using LTSS and Modified RSEI from 1990 to 2020—A Case Study of Zijin Mountain National Forest Park, Nanjing, China

Abstract

Evaluating the long-term urban forest ecological environmental quality (EEQ) and analyzing the drivers of its spatiotemporal changes can provide a scientific basis for making long-term urban forest planning decisions. Taking into account the characteristics of urban forest parks with low area proportions of construction land and bare land, high vegetation coverage, and serious forest disturbances, we constructed a modified urban forest park EEQ evaluation index based on a remote sensing ecological index named MRSEI, which is composed of the Landsat enhanced vegetation index (EVI), wetness, land surface temperature (LST), and forest disturbance index (FDI). We selected the Nanjing Zijin Mountain National Forest Park as the study area, used landsat time series stack (LTSS) remote sensing images from 1990 to 2020 as the main data source, and adopted the suggested modified MRSEI, the Theil-Sen median method, and the Hurst index to evaluate the EEQ to analyze its spatiotemporal variations and its driving factors in the study area. The main research results were as follows: (1) the EEQ of Zijin Mountain showed an up-and-down, overall slowly increasing trend from 1990 to 2020, while the spatial auto-correlation coefficient showed an overall decreasing trend; (2) the area percentage of the EEQ-persistent region accounted for 78.69%, and the anti-sustainable region accounted for 21.31%; (3) the spatial centers of the EEQ in the study area were mainly concentrated on the middle and upper part of the southern slope of Zijin Mountain, moving southward from 1990 to 2020; (4) the analysis of drivers showed that climate factors, forest landscape structure, forest disturbances, and forest growth conditions were the main driving factors affecting the EEQ in the study area. These results provide a research framework for the analysis of EEQ changes over a long-term period in the urban forest parks of China.

1. Introduction

The ecological environment is the basic guarantee for human survival and the foundation for social development. Dynamic changes in the ecological environment will directly affect the quality of human life [1]. Ecological environment quality (EEQ) generally refers to the quality of the ecological environment under the resources and environmental factors within a specific time and space range [2]. Studying the long-term regional EEQ changes and its driving factors can provide a scientific basis for the development of long-term regional forest planning.

With the advantages of efficient, real-time, and dynamic monitoring, remote sensing has become an effective method for assessing changes in regional forest ecosystems by providing a better description of ecological patterns and processes [3]. To a large extent, remote sensing can overcome these limitations. On 1 September 2008, the United States Geological Survey (USGS) opened Landsat archived remote sensing data to users worldwide free of charge. It is characterized by wide coverage, sufficient spectral information, high revisiting frequency, and spatial resolution (30 m× 30 m) similar to the size of forest management unit (1000 m²). Therefore, by using Landsat time series data (LTSS), we can provide a fine-scale description of regional forest ecosystem change and its driving factors on a regional long-time scale [4], which has become a hot-spot for remote sensing study in forestry [5].

At present, many methods including the analytical hierarchy process (AHP), fuzzy clustering method (FCM), and index evaluation method are widely used to evaluate the urban forests’ regional EEQ [6,7]. However, the above methods are somewhat subjective in the selection of indicators, making it difficult to comprehensively evaluate the ecological quality of the region [8]. The Ministry of Ecological Environment of China promulgated the “Technical Specification for Eco-environmental Status Evaluation” in 2006, which gives the definition and calculation method of the ecological index (EI) [9]. A total of five evaluation indices were selected in the specification, namely: biological abundance, vegetation coverage, water network density, land degradation, and environmental quality. The ecological environment index (EI) is formed through weighted summation of five indices. Since 2006, EI has been widely used in China, but a number of problems have been found in the trial, such as the reasonableness of the weights [10,11], the setting of the normalization coefficients [12,13], the ease of access to the indexes, etc. There is also an obvious shortcoming when the EI indexes can only illustrate the ecological status of an area in a general manner, but cannot be visualized to show the spatial distribution of different ecological conditions. In addition, spatial changes in the ecological environment cannot be analyzed. EI has coupled several ecological impact elements such as human activities, pollution indicators, biological abundance, water network density, vegetation cover, etc. However, the persistence of parameter weights makes it difficult to meet the monitoring needs of different ecosystems [14].

The greenness, humidity, heat, and dryness in the remote sensing ecological index (RSEI) proposed by Xu et al. are four important indicators in nature and closely related to human survival [15]. In RSEI, the greenness index is similar in meaning to the vegetation cover index in EI, and it is highly related to its biological abundance index; because the calculation bases of biological abundance and vegetation cover in EI are the same, only the weights are slightly different. The humidity indicator is the same as the water network index indicator in EI, but the humidity indicator can represent not only the open water body, but also the humidity of the soil and vegetation which is highly relevant to the ecological environment. The bare soil index indicator represents dryness and is closely related to the land degradation indicator in EI; the higher the bare soil index, the barer the surface, and the more serious the land degradation. Therefore, the three indicators in the new RSEI are closely related to the four indicators in the EI. Since the environmental quality index in the EI relies on annual statistical data rather than spatial data and is given the least weight in the EI, the index is not used in EI ecological grading. Therefore, the proposed RSEI index and ecological index (EI) are strongly comparable. Since the RSEI index established by principal component analysis (PCA) of the four indicator factors have a good agreement with the EI index, RSEI can be used as a new index for EEQ evaluation. RSEI has shortened the interval of EEQ monitoring with the advantages of fast calculation speed, less auxiliary parameters required, and high reliability of results.

Being located in urban planning areas, the EEQ of China’s urban forest parks is characterized by rapid spatial-temporal changes compared to that of suburban forests [16]. In the urban forest parks of southern China, the spectral signals of land conversion, forest cutting, forest tending, and regeneration on remote sensing images become unrecognizable after only a few years because of regional favorable hydrothermal conditions and short-term forest renewal cycles. At present, the ecological environment evaluation of forest parks based on RSEI in China usually takes long interval 1 phase or 2–3 phases of remote sensing data for static or dynamic evaluation [17]. It is difficult to capture the effects of burst forest disturbances (forest harvesting, deforestation), slow disturbances (climate change), and forest restoration processes (artificial regeneration, closing the mountains for reforestation) using this evaluation method. In addition, the studies on EEQ evaluation rarely analyze the long-term spatial-temporal changes and the driving factors of EEQs from the perspectives of forest disturbance, landscape pattern changes, and climate factors, which makes it difficult to delineate the complex relationships among forest disturbance, landscape pattern changes, climate change, and EEQs. Thus, these results cannot provide a scientific basis for the making of sustainable management measures for urban forest.

In China, as the green lungs of cities, urban forest parks usually have higher vegetation coverage [18]. At the same time, because most urban forest parks are located in urban planning areas, they are seriously affected by natural and human disturbances including tree cutting, forest land reversion, and the invasion of exotic pests [19]. Among the current EEQ evaluation studies based on the RESI, a set of indicator system targeting the characteristics of urban forest parks has not yet been formed [20]. In the RSEI model, the normalized difference building-soil index (NDBSI) used to express dryness consists of two parts: the building index and the bare soil index, which reflects the impact of urbanization on the regional land use pattern [21]. However, the calculation formula of NDBSI is complex, and the indicator is less targeted for urban forest parks with low building area and high vegetation coverage. In the RSEI model, greenness is usually expressed by the normalized difference vegetation index (NDVI), but the NDVI has the disadvantage of being easily saturated and having reduced sensitivity to areas with high vegetation coverage [22]. Compared with NDVI, the Landsat enhanced vegetation index (EVI) corrects for some atmospheric conditions and canopy background noise, and is more sensitive in densely vegetated areas of urban forest parks [23]. Compared with the NDBSI, the forest disturbance index (FDI) proposed by Masek et al. is relatively simple and can better distinguish forest pixels from other land use types [24]. Therefore, we constructed a modified RSEI called MRSEI by replacing NDVI with EVI and NDBSI with FDI.

Zijin Mountain National Forest Park is located in the eastern part of Xuanwu District in the Nanjing metropolitan area. Since the beginning of 1990s, as the park is located within the urban planning area, with convenient transportation and many residents living around, the ecological quality of the park has been disturbed by many natural and human factors such as invasion of alien pests and diseases, nibbled forest land by real estate development, environmental improvement projects around the periphery of the scenic area and the trampling of vegetation by mountaineering enthusiasts [25]. The ecological and environmental problems encountered in the Nanjing Zijin Mountain National Forest Park are very common in China, so they can be regarded a typical representative of urban forest parks in the country. In this study, we took the Nanjing Zijin Mountain National Forest Park as the case study area, and used the LTSS remote sensing images from 1990 to 2020 as the main information source, and adopted the MRSEI to analyze the spatial-temporal dynamics of the EEQ in the study area and its driving factors of spatiotemporal changes. Therefore, the objectives of this study were the following: (1) construct a modified RSEI which can specifically evaluate the EEQ of urban forest parks over 30 years in China; (2) reveal the spatial-temporal patterns of EEQ changes in the study area in long time series and predict the future trends of EEQ changes; (3) clarify the main driving factors affecting EEQ and provide suggestions on continuously improving EEQ of the urban forest park.

2. Materials and Methods

2.1. Study Area

Zijin Mountain National Forest Park, which is located in the eastern part of Nanjing, Jiangsu, China (Figure 1), is the first large urban forest park in China entitled “National Forest Park”. The park belongs to the first batch of national 5A-level tourist attractions, and is a demonstration site of Chinese ecological culture [26], playing a key role in the urban tourism of Nanjing. The park is divided into six divisions: the Mingxiao Mausoleum Scenic Division, the Sun Yat-sen Mausoleum Scenic Division, the Linggu Temple Scenic Division, the Mountain Park Scenic Division, the White Horse Park Scenic Division, and the Nature Landscape Scenic Division. In the park, Toutuo Ridge, with an altitude of 448.9 m, is the highest peak in the Nanjing and Zhenjiang mountain range of Jiangsu Province. Zijin Mountain is located on the edge of the northern subtropics, with annual precipitation of 900 to 1000 mm, and the average annual temperature of 15.7 °C [27].

The park has the unique advantages of superior natural conditions and abundant historical sites. There are more than 200 natural and historical sites in the area, which include one world heritage site, 15 major historical and cultural sites protected at the national level in China, and 31 cultural sites protected at the provincial level [28]. The forest cover rate in this study area is 78.2%, accounting for 15.6% of the total forest area in Nanjing [29]. With the development of the regional economy and the improvement of the living standards of local residents, Zijin Mountain National Forest Park has gradually changed from a suburban scenic area to a recreational urban park for the local residents. With the increasing spiritual and cultural needs of the tourists and local residents, the impact of various natural and human disturbances on the forest ecological environment has gradually increased. In recent years, the number of visitors to the forest park has been on the rise year by year, and the overall annual flow of tourists from other places has been around 10 million, while the annual number of local residents visiting the park fluctuating around 9 million [30]. In 1982, the alien invasive species Bursaphelenchus xylophilus was introduced to Zijin Mountain from Japan, and the area of forests infected by this alien pest has been expanding since then. A large number of scenic forest trees, which are mainly black pine and Masson pine, have died. At the same time, due to the impact of tourism activities, much of the wildlife, including Luehdorfia chinensis, a kind of rare butterfly which is on the list of National Grade II protected species, is at risk of extinction. The vegetation on both sides of the Zijin Mountain hiking trail are trampled and severely damaged. In 2009, a large area of landslides and serious soil erosion occurred on the southern slope of Zijin Mountain [31].

2.2. Data Collection and Pre-Processing

The data used in this study are: (1) Landsat TM/ETM+/OLI data for the study area from 1990 to 2020 provided by the USGS (http://earthexplorer.usgs.gov/, accessed on 1 July 2023). Image acquisition date is selected from April to October, which are the vegetation growing seasons. Cloud cover of most images is less than 5%. The spatial resolution of the multi-spectral band is 30 × 30 m and the panchromatic band is 15 × 15 m. Due to limitations of cloud cover and image striping, the images are selected in 2-year intervals (

Table 1. Basic characteristic of Landsat TM/ETM+/OLI data during 1990 to 2020.

Image ID	Acquisition Date	Sensor Type	Cloud Coverage/%	Data Level
LT51200381990304HAJ00	31 October 1990	Landsat TM	0.02	L1TP
LT51200381992294BJC01	20 October 1992	Landsat TM	11.86	L1T
LT51200381994139HAJ00	19 May 1994	Landsat TM	0.00	L1T
LT51200381996113HAJ00	12 April 1996	Landsat TM	0.00	L1T
LT51200381998150ULM00	30 May 1998	Landsat TM	0.00	L1T
LT51200382000284BJC00	10 October 2000	Landsat TM	0.00	L1T
LT51200382002193BJC00	12 July 2002	Landsat TM	0.00	L1T
LE71200382004127EDC02	6 May 2004	Landsat ETM+	0.00	L1TP
LT51200382006092BJC00	2 April 2006	Landsat TM	0.00	L1T
LE71200382008122EDC00	1 May 2008	Landsat ETM+	1.00	L1TP
LE71200382010095EDC00	5 April 2010	Landsat ETM+	1.00	L1TP
LE71200382012293EDC00	19 October 2012	Landsat ETM+	0.00	L1TP
LC81200382014162LGN01	11 June 2014	Landsat OLI_TIRS	8.86	L1TP
LE71200382016256EDC00	12 September 2016	Landsat ETM+	2.00	L1TP
LC81200382018301LGN00	28 October 2018	Landsat OLI_TIRS	0.06	L1TP
LE71200382020139EDC00	18 May 2020	Landsat ETM+	4.00	L1TP

); (2) “Detailed Planning of Sun Yat-sen Mausoleum Scenic Area” prepared by the Urban Planning and Design Research Institute of Southeast University, China; (3) “Peripheral Planning of Zijin Mountain Scenic Area” prepared by EDAW Landscape Design, San Francisco, United States in August 2004; (4) “Retrospective Environmental Impact Report on the Periphery of Sun Yat-sen Mausoleum Scenic Area” prepared by Nanjing Zijin Mountain Scenic Area Construction Development Co. Ltd, Nanjing, China; (5) “The Master Plan of Nanjing Zijin Mountain National Forest Park” prepared by Nanjing Forestry University (2014–2023).

Since the acquired remote sensing images of the study area from 1990 to 2020 belong to L1T/L1TP products and have been geometrically orthorectified by the USGS EROS Center before downloading, the pre-processing of Landsat TM/ETM+/OLI time series data is mainly focused on radiometric correction, water, and cloud detection [32,33]. The radiometric calibration tool in ENVI 5.3 was used to perform radiometric calibration of LTSS images. Then, the FLAASH tool of ENVI 5.7 software was applied to conduct the atmospheric correction of LTSS images. Automatic identification and masking of clouds was performed by using the algorithm for automatic cloud assessment [34]. There are many stripes on Landsat-7 ETM+ SLC OFF images acquired after 31 May 2003, so the add-on tool DeStripe [35] of ENVI was used to repair them.

After pre-processing, the LTSS was classified by using oriented classification software eCognition Developer 9.01 into six land use types: forest, crop, water, grass, built-up, and bare land for the temporal trajectory analysis of the forest landscape pattern index. A total of 350 random points were generated in ArcGIS based on the historical statistical area proportions of various land cover types in Purple Mountain National Forest Park in the “Nanjing Xuanwu District Statistical Yearbook”, including 150 forest points, 50 crop points, 40 water points, 30 grass points, 50 built-up points, and 30 bare land points. A total of 70% of the random points of each land type were randomly selected for training and 30% were used for testing. After supervised classification of the Landsat images, combined with the original remote sensing images, a visual interpretation was carried out with reference to the historical Google Earth images and field survey data of the corresponding year, and the actual land cover type of random point location was determined. Then, the overall classification accuracy and the classification accuracy of each land cover type were calculated accordingly. The calculation results show that from 1992 to 2020, the average overall classification accuracy of 16 remote sensing images was 89.38%, and the average classification accuracy of forest, crop, water, grass, built-up, and bare land was 91.84%, 83.40%, 96.00%, 80.76%, 99.45%, and 82.34%, respectively.

2.3. MRSEI

The modified remote sensing ecological index (MRSEI) is constructed using these four in dictators, as shown in Equation (1):

$$MRSEI-f\left(EVI,Wet,LST,FDI\right)$$

(1)

where EVI is the vegetation index, Wet is the soil moisture, LST is the surface temperature, and FDI is the forest disturbance index.

2.3.1. Landsat Enhanced Vegetation Index (EVI)

EVI is similar to NDVI and can be used to quantify vegetation greenness. However, EVI corrects for some atmospheric conditions and canopy background noise and is more sensitive in areas with dense vegetation [25]. As shown in Equation (2):

$$EVI=2.5\times \left(\frac{NIR-R}{NIR+6.0\times R-7.5\times B+1}\right)$$

(2)

where NIR/R/B is the near-infrared, red, and blue band of Landsat remote sensing images reflectivity.

2.3.2. Wetness

K-T transform is an empirical linear orthogonal transformation of images based on the information distribution structure of soil and vegetation in multidimensional spectral space of multispectral remote sensing, which is an effective data compression and redundancy removal method [36]. The brightness, greenness, and wetness components extracted from K-T transform are directly related to the physical parameters of the ground surface, so they have been widely used in ecological environment monitoring [37,38,39]. The humidity index used in this paper is represented by the wetness. However, the wetness equations vary with different sensors, as shown in Equation (3) by Chander et al. [40] for Landsat TM+ data; Equation (4) by Crist, Eric [41] for Landsat ETM+ data; and Equation (5) by Baig et al. [42] for Landsat OLI data:

$$We{t}_{TM+}=B\times 0.0315+G\times 0.2021+R\times 0.3102+N\times 0.1594-M_1\times 0.6808-M_2\times 0.6109$$

(3)

$$We{t}_{ETM+}=B\times 0.2626+G\times 0.2141+R\times 0.0926+N\times 0.0656-M_1\times 0.7629-M_2\times 0.5388$$

(4)

$$We{t}_{OLI+}=B\times 0.1511+G\times 0.1973+R\times 0.3283+N\times 0.3407-M_1\times 0.7117-M_2\times 0.4599$$

(5)

where B, R, G, M₁, and M₂ represent the blue, green, red, near-infrared band, mid-infrared 1, and mid-infrared 2 reflectance of TM/ETM+/OLI data, respectively.

2.3.3. Land Surface Temperature (LST)

To reduce the influence of different acquisition dates of remote sensing images on the study results, the remote sensing image on 2 July 2002 was chosen as the reference image, and 37 pseudo invariant feature (PIF) points from the buildings and roads in the study area were generated to perform relative corrections on other 15 phases of land surface temperature raster extracted from thermal infrared band, respectively [43]. Land surface temperature (LST) was calculated by using the model in Landsat user manual [44]. As shown in Equation (6):

$$T=\frac{K_2}{\ln \left(K_1+1\right)}$$

(6)

where T is the temperature of the sensor; K₁ and K₂ are two calibration coefficients. For TM+, K₁ = 60.776 mWcm⁻² sr⁻¹μm⁻¹, K₂ = 1260.56 K. For ETM+, K₁ = 66.609 mWcm⁻² sr⁻¹μm⁻¹, K₂ = 1282.71 K. For TIRS, K₁ = 77.489 mWcm⁻² sr⁻¹μm⁻¹, K₂ = 1321.08 K.

Sensor temperature T is converted to land surface temperature LST by using Equation (7) [45].

$$LST=T/\left[1+\left(\lambda T+\rho \right)\ln \epsilon \right]$$

(7)

where LST is the land surface temperature; λ is the center wavelength of the thermal infrared band; ρ = 1.438 × 10⁻² mK; ε is the emissivity of ground objects, and its value is estimated using NDVI based on Sobrino’s model [46].

2.3.4. Forest Disturbance Index (FDI)

In the histogram of the red band of Landsat image after masking dark matters, pixels whose reflection values are lower than the forest peak value are called pure forest pixels. On the basis of extracting LTSS pure forest pixels, tasseled cap transform was conducted for each image [26]. Forest disturbance index is calculated by Equation (8):

$$FD{I}_p=\frac{b_p-\overline{b}}{S{D}_p}-\frac{g_p-\overline{g}}{S{D}_g}-\frac{w_p-\overline{w}}{S{D}_w}$$

(8)

where $S{D}_p$, $S{D}_g$, and $S{D}_w$, respectively, represent the mean and standard deviation of brightness, greenness, and wetness of pure forest training samples on remote sensing images of the study area, while b_p, g_p, and w_p represent the brightness, greenness, and wetness of pixel p on the image. Since forest pixels have lower brightness values, higher greenness, and wetness values, it can be seen from Equation (8) that the smaller the FDI_p, the more likely the pixel is a forest pixel. Conversely, the larger the FDI_p, the more likely the pixel is a non-forest pixel.

2.3.5. Construction of MRSEI Model

Due to the dimensional inconsistency of the four indicators, the weights of each indicator would be unbalanced if used directly in the principal component analysis (PCA). It was necessary to standardize these four indicators before doing the analysis of PCA to unify their magnitudes between [0, 1]. Then, PCA was performed to obtain the variance contribution of each principal component. In PCA, most of the normalized indicators are explained by the first principal component (PC1) and then the initial remote sensing ecological index (MRSEI₀) was obtained [47]. The equation is shown in Equation (9):

$$MRSE{I}_0=1-\left\{PC1\left[f\left(EVI,Wet,LST,FDI\right)\right]\right\}$$

(9)

Additionally, the MRSEI₀ was standardized to facilitate index comparison. The closer the MRSEI₀ is to 1, the greater the EEQ [19,22].

2.3.6. Water Masking

Before the analysis of PCA, the water bodies needed to be masked because of the influence of the water bodies on the humidity indicator. The newly modified normalized difference water index MNDWI, which can solve the difficult problem of separating water and shadow, was used in this study. MNDWI can reduce the error of expanding the actual water body extent by using NDWI, and can better reflect the detailed characteristics of water, such as the content of suspended sediment [48]. MNDWI is calculated by Equation (10):

$$MNDWI=\frac{Green-MIR}{Green+MIR}$$

(10)

where MIR is the mid-infrared band, such as band 5 of the TM/ETM+ image, and Green is the green band, such as band 2 of the TM/ETM + image.

2.4. Analysis of Driving Factors

2.4.1. Landscape Pattern Index

Landscape pattern index can reflect the structure composition and spatial configuration of the research object at the landscape scale. Since highly condensed landscape structure information can be extracted from landscape pattern index, it is possible to measure the correlation between ecological process and spatial structure using the landscape pattern index [49]. The selection of a reasonable landscape pattern index has an important impact on forest landscape pattern research. Based on LTSS supervised classification, 72 forest patch indices from the aspects of average forest patch area, patch size, patch shape, patch connectivity, patch similarity, and core patch ratio were selected.

2.4.2. Analysis Method of Driving Factors

In this study, the EEQ extracted from MRESI₀ was selected as the dependent variable, and the mean annual temperature, mean annual precipitation, temperature in the month of image acquisition (month temperature), precipitation in the month of image acquisition (month precipitation), 72 landscape pattern indices of forest patches and vegetation index of greenness extracted from K-T transform were used as independent variables for Pearson correlation analysis. In addition, in order to analyze the influence of different factors on the EEQ, Pearson correlation analysis was used to analyze the driving factors. From the Pearson correlation analysis, the variables with significant correlation with the dependent variable were selected to establish the multiple regression model [50].

2.5. Analysis of Spatial-Temporal Variation of EEQ

2.5.1. Analysis of Temporal Variation

The Theil-Sen median, which is a calculation method applicable to trend analysis for long-term series data, is computationally efficient, insensitive to measurement errors. This method can be combined with the Mann–Kendall trend analysis method to study the trend of long-term EEQ change of Zijin Mountain [51,52]. The method is shown in Equation (11):

$$\beta =median\left(\frac{MRSE{I}_i-MRSE{I}_j}{i-j}\right)$$

(11)

where median represents taking the median value; MRSEI_i and MRSEI_j are the corresponding data of MRSEI in year i and j, respectively, and j > i; when β = 0, MRSEI has no change, when β > 0, MRSEI increase, and vice versa. The significance of the β values was further tested using the Mann–Kendall method [53].

2.5.2. Spatial Variation Analysis

Spatial auto-correlation analysis can reveal the regional structure pattern of spatial variables and obtain the spatial distribution characteristics. Global spatial auto-correlation describes the spatial dependence of an object within the total spatial range, which can be expressed by Moran’s I index and Geary‘s C index. Among these two indices, Moran’s I is more often used to verify whether adjacent areas within the study area have the same change trend (positive spatial correlation), opposite (negative spatial correlation), or independent (random distribution) [54,55]. Moran’s I > 0 means positive spatial correlation, and the larger the value, the more obvious the spatial correlation; Moran’s I < 0 means negative spatial correlation, and the smaller the value, the greater the spatial difference; otherwise, Moran’s I = 0 means random distribution.

The shift in the centers of gravity of EEQ can reflect the spatial and temporal aggregation of EEQ in the process of spatial evolution. The gravity center coordinates of weighted average EEQ (that is, the average of the centroid coordinates of all EEQ sample elements) are calculated. Then, the iterative EEQ is used to find the coordinate position of the point which can make the Euclidean distance between all EEQ in the data set reach the minimum. The gravity centers equation of EEQ is expressed as follows:

$$\overline{X_t}=\frac{{{\displaystyle \sum }}_{i=1}^n{w}_i{x}_i}{{{\displaystyle \sum }}_{i=1}^n{w}_i}, \overline{Y_t}=\frac{{{\displaystyle \sum }}_{i=1}^n{w}_i{y}_i}{{{\displaystyle \sum }}_{i=1}^n{w}_i}$$

(12)

where $\left(\overline{X_t},{\overline{Y}}_t\right)$ is the weighted average center coordinate; $\left(x_i,y_i\right)$ is the coordinate of EEQ sample elements; n is the sum of EEQ sample elements; and w_i stands for the weight of EEQ of sample i.

2.6. Future Change Trend Analysis of EEQ

The Hurst index based on the re-scaled range (R/S) analysis is an important indicator for quantitatively describing the long-term dependence of time series data, and its basic principle, is that within a time series ξ(t) (t = 1, 2, …), for any positive integer ≥ 1, there is a mean series [56]:

$${\zeta }_{\tau }=\frac{1}{\tau }{\displaystyle \sum }_{t=1}^{\tau }\xi \left(t\right) \left(\tau =1,2,\dots \right)$$

(13)

The cumulative deviation is obtained as:

$$X\left(t,\tau \right)={\displaystyle \sum }_{t=1}^{\tau }\left[\zeta \left(u\right)-{\zeta }_{\tau }\right]? \left(1\le t\le \tau \right)$$

(14)

The polar difference R is defined as:

$$R\left(\tau \right)=\underset{1\le t\le \tau }{\max }X\left(t,\tau \right)-\underset{1\le t\le \tau }{\min }X\left(t,\tau \right)? \left(\tau =1,2,\cdots \right)$$

(15)

The standard deviation S is defined as:

$$S\left(\tau \right)={\left[\frac{1}{\tau }{\displaystyle \sum }_{t=1}^{\tau }{\left(\xi \left(t\right)-{\xi }_{\tau }\right)}^2\right]}^{\frac{1}{2}} \left(\tau =1,2,\dots \right)$$

(16)

R, S, and τ satisfy the general relationship:

$$\frac{R\left(\tau \right)}{S\left(\tau \right)}=c\cdot {\tau }^H$$

(17)

Least squares fitting is defined as:

$$\log {\left(R/S\right)}_{\tau }=\log c+H\log \tau$$

(18)

where c is a constant, R(τ)/S(τ) is the re-scaled range, and H is the Hurst index. The Hurst index indicates the range of the cumulative deviation from the mean of the time series over time. The value of Hurst index of the time series is between 0 and 1. There are three main scenarios: (1) 0 < H < 0.5, indicating that the future trend of the EEQ is anti-persistence, i.e., it opposite to the past trend; (2) 0.5 < H < 1, indicating that of the future trend of the EEQ is persistence, i.e., it consistent with the past trend; or (3) H = 0.5, indicating that the future trend of the EEQ is random.

3. Results

3.1. Change Analysis of Annual Mean EEQ

In this study, the principal component transformation (PCA) was used to synthesize MRSEI which is represented by EVI, Wet, LST, and FDI from 1990 to 2020 with the information of the first principal component (PC1). Since PC1 had the largest covariance eigenvector among the four principal components and its contribution was above 75% (

Table 2. Principal component load matrix table from 1990 to 2020.

Year	Variable	PC1	PC2	PC3	PC4
1990	Covariance eigenvalue (%)	73.54	17.23	7.84	1.39
1992	Covariance eigenvalue (%)	79.81	14.53	4.76	0.91
1994	Covariance eigenvalue (%)	82.13	13.89	2.88	1.11
1996	Covariance eigenvalue (%)	80.27	11.13	7.04	1.55
1998	Covariance eigenvalue (%)	86.49	10.78	2.28	0.45
2000	Covariance eigenvalue (%)	85.52	11.62	2.31	0.55
2002	Covariance eigenvalue (%)	86.91	10.20	2.30	0.60
2004	Covariance eigenvalue (%)	82.47	12.48	4.33	0.71
2006	Covariance eigenvalue (%)	81.98	13.15	3.96	0.92
2008	Covariance eigenvalue (%)	82.66	10.89	5.78	0.67
2010	Covariance eigenvalue (%)	77.32	12.31	8.60	1.78
2012	Covariance eigenvalue (%)	75.52	15.92	7.47	1.09
2014	Covariance eigenvalue (%)	82.65	13.87	3.16	0.32
2016	Covariance eigenvalue (%)	77.83	14.91	6.51	0.75
2018	Covariance eigenvalue (%)	76.34	18.56	4.75	0.34
2020	Covariance eigenvalue (%)	78.53	15.76	5.21	0.50

), it was figured out to represent the EEQ in this study according to Equation (9).

The closer the MRSEI value is to 1, the better the EEQ. The MRSEI values were divided into five equal-interval classes (Figure 2): “excellent” (0.8–1.0); “good” (0.6–0.8); “medium” (0.4–0.6); “poor” (0.2–0.4); and “bad” (0.0–0.2) [18,21]. The results showed that from 1990 to 2020, the overall MRSEI of Zijin Mountain was high and the overall EEQ was good, while the small areas with low EEQ were located in the southern and eastern part of the park with dense residents and convenient road traffic conditions.

From 1990 to 2020, the EEQ in the study area showed an overall slow upward trend, rising from 0.754 in 1990 to 0.762 in 2020. To further investigate the pattern of change, we divided the EEQ variations in the study period into 5 phases (Figure 3): the EEQ showed a steady upward trend from 1990 to 2002 (phase 1), a downward trend from 2002 to 2010 (phase 2), followed by an upward trend from 2010 to 2014 (phase 3), then a downward trend from 2014 to 2018 (phase 4), and finally a slow upward trend from 2018 to 2020 (phase 5). In 1990, with the deepening of reform and opening up, regional economic development was relatively rapid, and the Nanjing municipal government increased the investment in ecological environment construction. After the national market economy status was established in 1992, a real estate boom began to occur around the study area, and some forest land was converted into residential areas. Since 1993, a real estate boom has emerged across the country. To promote the tourism economy in Nanjing, many hotels, restaurants, and other tourist facilities such as an aquarium named the “Underwater World” were constructed in the park, causing the decrease in forest land, which affected the EEQ in the study area during 1994–1996. From 2002 to 2010 (phase 2), 13 villages and several factories on the peripheral areas of the park were moved out of the study area and replaced by many trees and grasses. During the process of the large-scale environmental improvement project named Relocation and Greening, the area of bare land caused by a large number of infrastructural projects increased dramatically, causing the decrease in EEQ in the park. During phase 3 (2010–2014), after the environmental improvement project was completed in 2008, the forest cover in the park increased, together with an increase in EEQ. From 2014 to 2018 (phase 4), in order to enrich the tourist resources, a number of tourism communities and cultural museums had been established, causing the increase in built-up area and the decrease in EEQ in the park. Finally, due to the outbreak of the coronavirus disease in 2019 (COVID-19), the number of tourists had decreased sharply, and human disturbance had also been greatly reduced, making the EEQ stably increase during 2018–2020 (phase 5).

3.2. Analysis of EEQ Spatial-Temporal Variation

3.2.1. Temporal Variation Analysis

The results of the Theil-Sen median slope estimation for the EEQ during the period of 1990–2020 indicated that the area of positive β values was significantly larger than the area of negative values (Figure 4). Specifically, there were 92.92% of the regions with an increasing trend and 7.08% of the regions with a decreasing trend. The mean value of β was 2.02, indicating that the EEQ of Zijin Mountain had been in an overall increasing trend in the past 31 years. We further performed Mann–Kendall test on the EEQ trends, and the results showed that there were significant spatial heterogeneity of the EEQ of Zijin Mountain (Figure 5). The area percentage with an extremely significant increase accounted for 39.21%, the non-significant increases accounted for 54.81%; the area with an extremely significant decrease accounted for only 0.72%, the non-significant decrease regions accounted for 5.26%, and the percentage of regions with no change was small. In terms of spatial distribution, the areas of significant increase were mainly located in the southwestern part of the region with gentle slope and many cultural attractions; the area with significant decrease was mainly distributed in the southeast part and the edge of the park with dense residents and convenient traffic conditions. Regions of non-significant decrease were intertwined with regions of significant decrease, mainly distributed in residential areas on the edge of the park.

3.2.2. Spatial Auto-Correlation Analysis

The values of the Moran’s I for EEQ during the period of 1990–2020 were given in Figure 6. As shown in Figure 6, the spatial auto-correlation coefficients of the EEQ in the park were all positive, with a mean value of 0.6532, indicating that the EEQ of Zijin Mountain had been in a distribution pattern of significant spatial aggregation from 1990 to 2020. The trend line showed an overall decreasing aggregation trend of the EEQ in the study area from 1990 to 2020. This trend was closely related to the continuous intensive forest management measures and greening projects on the edge of the park during the study period. During the period 2002 to 2008, a large-scale relocation and greening project on the edge of the park was launched. Many villages and factories were moved out of the park and a lot of trees and grasses were planted on the original building sites. During the process of relocation and greening, the area of bare land increased, causing the widening EEQ gap between the edge area and the core area of the park. In 2008, the project was finished and a large amount of construction land was replaced by forest and grass, rapidly reducing the EEQ gap between the edge and the core area.

3.2.3. Spatial Center Dynamic Analysis

The EEQ spatial center dynamic analysis was performed; the standard deviation ellipse and the 16-phase EEQ geographic gravity centers together with an ellipse were generated (Figure 7). It can be seen from Figure 7 that the geographic gravity centers of the EEQ of Zijin Mountain were mainly concentrated on the upper part of the southern slope, gradually shifting southward. The reason for the spatial distribution patterns of the EEQ centers was closely related to the distribution pattern of core scenic sites and large-scale project of relocation and greening on the edge of the park. The most famous scenic sites, such as the Sun Yat-sen Mausoleum, Mingxiao Mausoleum, and Linggu Temple, are all concentrated on the upper part of the southern slope. Generally, the forest cover and quality around these famous scenic sites are higher than other places. Since the beginning of relocation and greening project on the southern edge of the park, a large area of built-up land had been replaced by trees and grasses, resulting in the moving of the geographic gravity center of the EEQ from the upper part to the middle and lower part of the southern slope. The center points of 2004 and 2010 are outside the ellipse, which may be related to the fallowing of farmland in 2004 and the beautification and greening of the environment in 2010.

The major semi-axis of the ellipse represents the direction of data distribution, and the minor semi-axis represents the range of data distribution. The greater the value difference between the long and short semi-axes (the greater the oblateness), the more obvious the directionality of the data.

3.3. Future Change Trends of EEQ

The future change trend of the EEQ in the park was analyzed at the pixel level using the Hurst index from 1990 to 2020 (Figure 8). The Hurst index was divided into five classes according to the value range: strong anti-persistence (SAP) (≤0.25), weak anti-persistence (WAP) (0.25 < H ≤ 0.5), random (R) (H = 0.5), weak persistence (WP) (0.5 < H ≤ 0.75), and strong persistence (SP) (H > 0.75). The area percentage of the EEQ persistent region in the park was 78.69%, and the anti-persistent region was 21.31%. Among the four classes, the weak persistence ratio was the biggest with a value of 32.26%, revealing that the future EEQ in the study area would continue to improve, i.e., with the same trend as the past.

In order to explore the persistence of the change trend of the EEQ in the park, the future change trend of the EEQ was obtained by combining the Hurst index and the grading results of the trend test. A total of 8 trends were produced (Figure 9). Among them, the percentages of non-significantly increasing-weakly anti-persistent (NSIWAP), non-significantly increasing-weakly persistent (NSIWP), significantly increasing-weakly anti-persistent (SIWAP), and significantly increasing-weakly persistent (SIWP) were 28.53%, 32.85%, 15.23%, and 5.29%, respectively. The large area percentage of non-significantly increasing-weakly persistent (NSIWP) indicated that the EEQ of the study area may exhibit an increasing trend in the future. This trend is due to the long-term intensive management measures of the scenic forest such as the relocation and greening project on the edge of the park, which have played a very important role in the improvement of the EEQ of the study area.

3.4. Analysis of EEQ Driving Factors

Studies have shown that EEQ in forest areas is related to vegetation growth condition, forest disturbance, and climatic factors [57]. At the landscape scale, EEQ is also closely related to forest landscape structure. Therefore, the EEQ during 1990 to 2020 was used as the dependent variable, the greenness extracted from the K-T transform was used as an index reflecting the vegetation growth condition, 72 landscape pattern indicators, and four climate factors of annual average temperature, annual average precipitation, the precipitation in the month of image acquisition (month precipitation), and the temperature in the month of image acquisition (month temperature), were used as independent variables to conduct Person correlation analysis. It can be seen from (

Table 3. Impact factors significantly correlated with EEQ.

Impact Factor	Correlation
Greenness	0.611 **
Month temperature	0.606 **
Patch similarity	−0.439 *
Patch connectivity	0.440 *
Forest Disturbance Index	−0.458 *

Note: * indicates significant at the 0.05 level, ** indicates significant at the 0.01 level.

) that greenness and month temperature are significantly correlated with EEQ (α = 0.01), and patch connectivity, patch similarity, and FDI are significantly correlated with EEQ (α = 0.05).

4. Discussion

This study constructs a modified RSEI, which can more accurately evaluate the EEQ of an urban forest park. During the calculation of PCA of MRESI, in order to avoid the influence of input sequence of indicator and number of principal components on EEQ, the four indicators of 16 phases were all input strictly according to the sequence of EVI, greenness, LST, and FDI. To reduce the influence of inconsistent image acquisition time on the results of LST, the land surface temperature on 12 July 2002, which was the highest during 1990 to 2020, was chosen as the reference to relatively correct LST of other periods using the pseudo invariant feature (PIF) method. In order to reduce the errors caused by different dimensions, the four indicators were normalized before the analysis of PCA. Since PC1 and PC2 also contain usable information, contribution-weighted linear combinations of the three principal components, PC1, PC2, and PC3, should be considered in future studies to improve the quality of EEQ analysis.

The research results indicate that from 1990 to 2020, the EEQ of the park shows an up-and-down, but slowly increasing trend, while the spatial auto-correlation coefficient of the EEQ shows an overall decreasing trend, and the future EEQ will continue to improve. The research results of scholars on the ecological environment quality of Nanjing city indicate that the EEQ in the Zijin Mountain area is gradually improving, which is consistent with the research results of this study [58]. The study by Zhang Peng et al. shows that the ecological environment quality of forest parks in Nanjing has changed dramatically in the past thirty years, and the trend of change is consistent with the results of our study [59]. These trends of change in the study area are closely related to the long-term intensive forest management measures from 1992 to 2020, including the relocation and greening project on the edge of the park [23]. Being in a national 5A scenic area and the leading role in municipal tourism industry, the forest resources in the park have been under strict protection for a long time, except for the edges with convenient transportation and high population density. Protected under such a situation, the forest ecosystems have been in a positive succession trend from coniferous forest to mixed forest of coniferous and broad-leaved, then to mixed forest of broad-leaved, resulting in sustained improvement of forest quality [60]. In order to improve the EEQ of the park, a large-scale relocation and greening project had been carried out on the edge of the scenic area since 2002, making the EEQ gap between the edge of the park and the internal core area gradually narrowed. On 22 September 2020, Chinese government officially proposed the goal of achieving carbon peaking by 2030 and carbon neutrality by 2060 at the 75th session of the United Nations General Assembly [61]. As the green lung of Nanjing City, the forest in the park plays an extremely important role in carbon sequestration and oxygen release, and has been paid unprecedentedly attention by local government and residents. Therefore, the forest management intensity will be further improved and the future EEQ will show a continuously increasing trend.

Forest vegetation conditions, climate factors, and the forest landscape structure are the main drivers affecting the long-term EEQ of the study area [62,63]. Since the forest cover of the study area is as high as 78.2%, the EEQ is closely related to the conditions of the forest growth conditions [24]. Temperature, as an important factor, is a key environmental factor affecting vegetation growth [64,65]. Studies have shown that EVI usually increases with the increase in temperature. Forest pests and diseases, conversion of forest land, afforestation projects, and other forest disturbances affect the quality of the regional EEQ by changing the forest area and affecting forest growth [23]. Changes in forest landscape patterns such as connectivity and similarity of forest patches often trigger changes in water–heat interactions at the “ground-atmosphere” interface in the study area, affecting microclimate features and thus changing forest growth conditions. It should be noted that the above three kinds of factors do not affect the EEQ of the study area in isolation, but are closely related to each other [11]. Changes in climate in the study area affect the growth of forest trees, further causing changes in forest landscape patterns, and ultimately affecting EEQ [66].

After China’s reform and opening up in 1978, with the acceleration of urbanization, many distant suburban parks became urban forest parks, and forest tourism relying on urban forest parks has become a pillar industry for regional economic development in many cities. With the improvement of traffic conditions, and the increase in tourists, many forest disturbances causing negative impact on the park have happened. These disturbances mainly include encroachment of forest land by real estate development, forest landscape fragmentation caused by tourist roads, invasion of alien forest pests and diseases, and trampling of forest vegetation by local mountaineering enthusiasts [19]. However, forest disturbance has both advantages and disadvantages. While reducing the EEQ of the park, active forest management measures such relocation and greening projects greatly improved the EEQ of fringe areas of the park, narrowing the EEQ gap between the edge area and internal core area. The ecological and environmental problems encountered in Nanjing Zijin Mountain National Forest Park are very common in China, especially in southern part of the country. The EEQ analysis method in urban forest parks proposed in this paper has the advantages of low cost, fast calculation, less auxiliary parameters required, and high reliability of results, and can be extended to other urban forest parks in southern China.

5. Conclusions

From 1990 to 2020, the EEQ of Zijin Mountain showed an up-and-down and overall-slowly increasing trend. Forest growth conditions, climate factors, forest landscape structures, and forest disturbances are the main driving factors affecting the long-term EEQ in the study area. Therefore, preventing and controlling invasive pests and diseases, reducing negative forest disturbances such as the conversion of forest land, improving forest growth conditions through forest tendering, and improving forest landscape quality by increasing forest patch connectivity are the main ways to improve the ecological quality in the study area.