Study on Regional Eco-Environmental Quality Evaluation Considering Land Surface and Season Differences: A Case Study of Zhaotong City

Abstract

Timely and quantitatively evaluating regional eco-environmental quality (EEQ) is of great significance for realizing regional sustainable development goals. Especially for cloudy areas, it was a great challenge to construct a regional EEQ dataset with high quality and high resolution. However, existing studies failed to consider the influence of land surface and season elements in evaluating regional EEQ. Therefore, this study aimed to promote an accurate EEQ-evaluating framework for cloudy areas. Zhaotong city, a typical karst and cloudy region, was chosen as the study area. First, we integrated multi-source spatiotemporal datasets and constructed a novel eco-environmental comprehensive evaluation index (ECEI) to assess its EEQ from 2000 to 2020. Next, standard deviation ellipse (SDE) and trend analysis methods were applied to investigate regional EEQ’s change trends. Finally, ecological index (EI) values for different years were calculated to validate the effectivity of the ECEI. The main findings were as follows: (1) The EEQ of Zhaotong showed an upward-fluctuating trend (0.0058 a⁻¹), with average ECEI values of 0.729, 0.693, 0.722, 0.749, and 0.730. (2) The spatial distribution pattern of the EEQ showed high values in the north and low values in the south, with Zhaoyang district having the lowest ECEI value. (3) From 2000 to 2020, the standard deviation of the major axis of the ellipse moved northeast of Zhaotong city with θ of SDE changing from 57.06° to 62.90°, thus, indicating the improvement of northeastern regions’ EEQ. (4) The coefficients of the determinant (R²) between the EI and ECEI were 0.84, which was higher than that of EI–RSEI (R² = 0.56). This indicated that our promoted framework and the ECEI could acquire more accurate EEQ results and provide suggestions for relevant policymakers.

1. Introduction

Since the Industrial Revolution, the urbanization of the whole world has increased dramatically [1]. Based on the United Nations report issued in 2014, the world’s urbanization rate is expected to reach 66% in 2050 [2]. Especially in China, since the implementation of the reform and opening-up policy, the urbanization rate has increased from 17.92% in 1978 to 64.72% in 2021, with an average annual growth rate of 103.03% [3]. However, with the rapid development of urbanization, a series of problems have arisen, for example, degradation of land quality and soil erosion [4,5], decrease in biodiversity and vegetation coverage [6,7], intensification of heat island effect [8,9], a decline of air quality [10,11], and deterioration of eco-environment [12]. In response to these problems, the Chinese government carried out a regional, coordinated development strategy in 2017 [13], aiming to decrease the gap between different regions and promote coordinated development between the eco-environment and urbanization. Specifically, it is necessary to improve these eco-environment issues from the perspective of green infrastructure and urban sprawl. Urban green infrastructure is defined as “a strategically planned network of natural and semi-natural areas with other environmental features designed and managed to deliver a wide range of ecosystems services” [14]. Marando et al. used European functional urban areas to study the influence of the implementation of urban green infrastructure in mitigating the urban heat island effect, and they found that a tree cover of at least 16% was required to achieve a l ℃ drop in urban temperatures [14]. Donati et al. found that the planning of blue-green infrastructure could promote the permeability and availability of “stepping stone” habitats in densely populated landscapes, indicating that the blue and green infrastructure construction was of great help in maintaining regional habitat connectivity [15]. When Chanchipricha and Fischer evaluated the impact of developing urban green infrastructure, they suggested that environmental impact assessment should consider green infrastructure [16]. Barbosa et al. concluded that there was a disconnect and disarticulation between public spaces and green areas [17]. In general, green infrastructure construction could provide some influence on promoting an urban eco-environment, including the mitigation of the urban heat island effect and the maintenance of regional habitat connectivity.

As for urban sprawl, its influence on the regional EEQ is complicated. Ewing and Hamidi found that sprawling areas had higher ozone levels than compact areas [18]. Polidoro et al. found that urban sprawl could lead to the consequences of territorial ordering and the egalitarian spatial distribution of essential services [19]. Tian et al. concluded that planning was strongly correlated with urban sprawl, which was kind of a “planned sprawl” [20]. Huang et al. took eight landscape metrics to analyze the spatiotemporal changes in landscape patterns, indicating that a guiding plan was required to integrate planning in both urban and non-urban planned districts [21]. Barbosa and Pradilla concluded that the degree of vulnerability encountered by the social urban spatial structure was higher in expansion areas than in central areas [22]. Generally, urban sprawl, to some extent, could influence the regional EEQ. Thus, it is important to integrate planning to guide the direction of urban sprawl. Therefore, to realize this goal, it is urgent to evaluate the regional eco-environmental quality (EEQ) timely and quantitatively.

The definition of an eco-environment is the total quantity and quality of water resources, land resources, biological resources, and climate resources that affect human survival and development [23]. Evaluation of the EEQ normally consisted of two parts, which were the construction of a scientific index system and an accurate assessment of the regional EEQ [24]. To date, many scholars have promoted numerous novel indexes to evaluate the regional EEQ. For example, Fu developed a comprehensive index system to evaluate the regional EEQ at the provincial scale; this index system included water and soil loss, land salinization, agricultural natural disasters, land resources, water resources, forest resources, air pollution, water pollution, solid waste pollution, population, economic development, and township enterprises. However, all used datasets were acquired from national statistical yearbooks [25]. He et al. combined fine particulate matter (PM_2.5) concentration, land surface temperature (LST), and vegetation cover (VC) to develop a comprehensive evaluation index (CEI), and he found that the increase of PM_2.5 concentration had a great influence on the regional EEQ. Specifically, degraded regions were mainly distributed in expanded urban areas with an increased PM_2.5 concentration [26]. Wei et al. integrated the normalized difference vegetation index (NDVI), wetness index, albedo, index-based built-up index, salinization index, and LST to construct an environmental quality index. In their study, weights of multiple indicators were objectively determined by the eigenvalues of all principal components. In addition, they found that the EEQ of urban regions displayed a damaging trend [27]. Chang et al. used the average annual precipitation, net primary productivity (NPP), elevation, biological abundance index (BAI), highway net density, and other nine indexes to establish an ecological environmental index (EEI). In their study, four climate aspects, such as soil–terrain, biological resources, and human factors were used to evaluate the regional EEQ. They also found that human factors had a negative influence on the regional EEQ [28]. Sun et al. calculated the NDVI, habitat quality index, normalized difference moisture index (NDMI), and normalized difference soil index (NDSI) to promote an eco-environmental quality index (EQI). Among them, threshold values were set to divide each index into different types. In addition, they also came to the conclusion that residential regions had poor EEQ [29]. Wei et al. selected twenty-three indices from hydro-meteorological, socioeconomic, soil, biological, and topographical aspects to construct an environmental vulnerability distance index (EVDI). Similar to the aforementioned studies, land use changes and human activities indirectly influenced the results of environmental vulnerability [30]. Chai and Lha chose sixty-six indexes from aspects of human activity and the natural environment to evaluate different regions’ EEQ; however, it took a long time to collect these data [31]. Ariken et al. built an index system from level, pollution, and protection aspects to evaluate the regional EEQ. Numerous indexes were calculated or acquired by statistical yearbooks [32]. To sum up, the data sources used in existing studies could be divided into two types: statistical data and non-statistical data [33]. Statistical data has the advantage of high authority. However, its data is mainly at the administrative scale, which fails to acquire the index value of any location [34]. Hence, in the previous studies, the regional EEQ was evaluated qualitatively. At the same time, it failed to detect EEQ changes at the grid level. Among non-statistical data, remote sensing data are one of the most widely used due to the advantages of wide space coverage and strong time continuity [35]. EEQ’s evaluation based on remote sensing datasets provided a useful way to detect EEQ changes at multi-dimensions, such as degradation of the EEQ caused by urban sprawl and human negative activities. To date, numerous indexes have been proposed to evaluate the regional EEQ, such as NDVI for vegetation detection [36] and LST for urban heat islands [37]. However, a single index could only reflect one aspect of the regional eco-environment, which fails to evaluate the regional comprehensive EEQ [38]. Compared with a single index, the aggregated index showed great potential. Especially, the remote sensing ecological index (RSEI) integrated four aspects, such as greenness, wetness, dryness, and heat, to evaluate the regional EEQ [39]. Since its promotion in 2013, numerous studies have applied this index to evaluate the regional EEQ at various scales [40,41,42,43,44]. However, the RSEI had some limitations in practical applications, including the influence of season, the effect of water bodies, lack of consideration for the land surface difference, and lack of application in cloudy areas. In response to the first three shortcomings, Xu et al. combined the land use data to develop the RSEI-2 to assess China’s EEQ [45]. Zhang et al. assessed Tianjin’s EEQ by combining the RSEI and season [46]. Xu and Deng suggested that it was not suitable to name these indexes as RSEI or use RSEI-related names when considering non-ecological factors [47]. Moreover, existing studies failed to consider both land surface and season differences. As for the fourth shortcoming, numerous studies adopted Google Earth Engine (GEE) platform to evaluate the regional EEQ by filtering satisfied images [48,49,50].

Generally, EEQ evaluation achieved fruitful results; however, cloudy regions’ EEQ evaluation considering both land surface and season differences has not been reported. This study aimed to promote an accurate EEQ-evaluating framework for cloudy areas and provide constant regional EEQ datasets with high quality and high resolution. Zhaotong city, a typical cloudy and karst area in China, was chosen as a study area to evaluate its EEQ from 2000 to 2020 as, during this period, the economic and social level of Zhaotong increased significantly. In addition, several high-accuracy datasets, such as land use datasets, have been available since 2000. First, Zhaotong’s cloudy-free Landsat series images from different seasons were filtered based on the GEE platform. Then, the results from the abundance index (AI) and four other indexes obtained in different seasons were integrated to establish a novel eco-environmental comprehensive evaluation index (ECEI). Finally, the ECEI results of were analyzed based on standard deviation ellipse (SDE) and trend analysis methods and compared with the ecological index (EI) developed by the Ministry of Ecology and Environment. This study is organized as follows: the Section 2 presents the general situations of the study area and adopted methodologies, which include calculation methods of the ECEI, EI, SDE, and trend analysis; the Section 3 describes the spatiotemporal distribution of the ECEI at different scales and the spatiotemporal change trend of the ECEI based on SDE and trend analysis methods; the Section 4 analyzes the spatial principal component analysis (SPCA) results, compares the ECEI results with the EI, and displays the implications and limitations of this study; the Section 5 shows this study’s main results and findings.

2. Materials and Methods

2.1. Study Area

Zhaotong is located in the northeastern part of Yunnan Province in China and the hinterland of the Wumeng Mountains on the border of Yunnan, Guizhou, and Sichuan provinces (Figure 1). It is also a typical karst region [51]. Zhaotong has a highland monsoon climate, with 96% of its territory covered by mountainous areas. It is rich in water resources, with an average annual precipitation of about 1120 mm and an uneven spatial and temporal distribution of rainfall, with the flood season concentrated in May–October. In 2020, the population and gross domestic product of Zhaotong were 5.09 million people and 128.87 billion yuan [52].

2.2. Data Sources and Processing

The study adopted multiple remote sensing datasets and administrative boundary data to evaluate Zhaotong’s EEQ. A brief introduction of each dataset is shown in

Table 1. Brief introduction of datasets.

Name	Resolution	Data Availability	Brief Description
Landsat series	30 m	https://www.usgs.gov/landsat-missions/landsat-collection-2-level-2-science-products accessed on 1 May 2022	A product of land surface spectral reflectance
GLC_FCS30	30 m	https://data.casearth.cn/ accessed on 15 May 2022	A product of global land cover with fine classification system
ASTER GDEM	30 m	http://www.gscloud.cn/home accessed on 22 May 2022	A product of global digital elevation model
Precipitation	1000 m	http://www.geodata.cn accessed on 5 July 2022	A product of monthly precipitation data
CHAP	1000 m	https://weijing-rs.github.io/product.html accessed on 20 June 2022	Products of China’s high-spatial air pollutants
HWSD	30 arc-second	https://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/ accessed on 20 June 2022	A product of harmonized world soil database
Statistical yearbook	\	http://www.zt.gov.cn/ accessed on 20 June 2022	Statistical data of regional socio-economic aspects
Administrative boundary data	\	http://www.ngcc.cn/ngcc/html/1/index.html accessed on 5 May 2022	A vector dataset for data mask and spatial analysis

.

Based on Table 1, Landsat 5, 7, and 8 images, land surface reflectance products of different seasons in the years 2000, 2005, 2010, 2015, and 2020 were obtained by using the GEE platform and median method. Due to the cloud pollution problem, it was difficult to get cloud-free images for each target year. Therefore, images from adjacent years were also included when filtering. Specifically, spring is from 1 March to 31 May, summer is from 1 June to 31 August, autumn is from 1 September to 30 November, and winter is from 1 December to 28 February [53]. GLC_FCS30 datasets were acquired from the Data Sharing and Service Portal website at a resolution of 30 m. ASTER GDEM data were obtained from the Geospatial Data Cloud website at a resolution of 30 m. Precipitation data were downloaded from the National Earth System Science Data Center with a 1000 m resolution. CHAP datasets were acquired from Weijing’s website at a resolution of 1000 m. HWSD data were derived from the Food and Agriculture Organization of the United Nations by selecting surface soil texture attributes with a 30-arc-second resolution. All raster datasets (Table 1) were resampled to 30 m resolution and converted to WGS_1984_48N coordinate system using ArcGIS 10.8 software.

2.3. Methods

The framework of this study is shown in Figure 2, which consists of three main parts: (1) ECEI construction; (2) EI construction; (3) ECEI accuracy assessment, spatiotemporal and trend analysis. Among them, SDE was used to present the change direction at the city scale. Section 2.3.1, Section 2.3.2, Section 2.3.3 and Section 2.3.4 describe the calculation formulas in detail.

2.3.1. Calculation of ECEI

The ECEI integrates NDVI, wetness (WET), normalized difference build-up and soil index (NDBSI), LST, and AI. After calculating each index for the four seasons, average ECEI values of different years were calculated. To reduce the influence of outliers, the values of the first 0.1% and last 0.1% in the histogram were removed. Equations (1)–(14) were utilized to obtain Zhaotong’s ECEI in 2000, 2005, 2010, 2015, and 2020, respectively [54,55]. All equations were calculated using the GEE platform and ArcGIS 10.8 software.

$$NDVI=\left({\rho }_{nir}-{\rho }_{red}\right)/\left({\rho }_{nir}+{\rho }_{red}\right)$$

(1)

$$\begin{array}{c}WE{T}_{TM}=0.03\times {\rho }_{blue}+0.20\times {\rho }_{green}+0.30\times {\rho }_{red}+0.16\times {\rho }_{nir}\\ -0.68\times {\rho }_{swir1}-0.61\times {\rho }_{swir2}\end{array}$$

(2)

$$\begin{array}{c}WE{T}_{ETM+}=0.26\times {\rho }_{blue}+0.21\times {\rho }_{green}+0.09\times {\rho }_{red}+0.07\times {\rho }_{nir}\\ -0.76\times {\rho }_{swir1}-0.54\times {\rho }_{swir2}\end{array}$$

(3)

$$\begin{array}{c}WE{T}_{OLI}=0.15\times {\rho }_{blue}+0.20\times {\rho }_{green}+0.33\times {\rho }_{red}+0.34\times {\rho }_{nir}\\ -0.71\times {\rho }_{swir1}-0.46\times {\rho }_{swir2}\end{array}$$

(4)

$$\begin{array}{c}NDBSI=\frac{1}{2}\left\{\frac{\frac{2\times {\rho }_{swir1}}{{\rho }_{swir1}+{\rho }_{nir}}-\frac{{\rho }_{nir}}{{\rho }_{nir}+{\rho }_{red}}-\frac{{\rho }_{green}}{{\rho }_{green}+{\rho }_{swir1}}}{\frac{2\times {\rho }_{swir1}}{{\rho }_{swir1}+{\rho }_{nir}}+\frac{{\rho }_{nir}}{{\rho }_{nir}+{\rho }_{red}}+\frac{{\rho }_{green}}{{\rho }_{green}+{\rho }_{swir1}}}\right\}+\\ \frac{1}{2}\left\{\frac{{\rho }_{swir1}+{\rho }_{red}-{\rho }_{nir}-{\rho }_{blue}}{{\rho }_{swir1}+{\rho }_{red}+{\rho }_{nir}+{\rho }_{blue}}\right\}\end{array}$$

(5)

$$LST=\left\{a\left(1-C-D\right)+\left[b\left(1-C-D\right)+C+D\right]\times T_b+D\times T_a\right\}/C$$

(6)

$$C=\epsilon \tau$$

(7)

$$D=\left(1-\epsilon \right)\left[1+\left(1-\epsilon \right)\tau \right]$$

(8)

$$T_a=16.01+0.09\times T_0$$

(9)

$$AI=\mu \times \left(\begin{array}{c}0.35\times Forest+0.21\times Grassland+0.28\times Water+\\ 0.11\times Cropland+0.04\times Build+0.01\times Unused\end{array}\right)/Area$$

(10)

$$X_{rescale}=\left(X_i-X_{\min }\right)/\left(X_{\max }-X_{\min }\right)$$

(11)

$$ECE{I}_{season\_o}=PC1\left\{f\left(NDVI,WET,NDBSI,LST,AI\right)\right\}$$

(12)

$$ECE{I}_{season}=\frac{ECE{I}_{season\_o\_i}-ECE{I}_{season\_o\_\min }}{ECE{I}_{season\_o\_\max }-ECE{I}_{season\_o\_\min }}$$

(13)

$$ECEI=\frac{{\displaystyle \sum _{j=1}^{N}ECE{I}_{season\_j}}}{N}$$

(14)

where ρ_blue, ρ_green, ρ_red, ρ_nir, ρ_swir1 and ρ_swir2 represent the land surface reflectance of the blue, green, red, nir, swir1, and swir2 bands, respectively; WET is the moisture component of the tasseled cap transformation, its calculation formulas are different for different Landsat series [56,57]; a and b are constant variables, which are −67.36 and 0.46, respectively; ε is the land surface emissivity; τ is the atmospheric transmittance; T_a is the average atmospheric action temperature; T_b is the brightness temperature, which is derived from Landsat images metadata file; T₀ is the actual surface temperature; μ is the normalized coefficient; Forest, grassland, water, cropland, built and unused are the corresponding areas of each land use type in the target units of the study area; X_rescale indicates the normalized result of each index; X_i, X_min and X_max represent the i_th, min and max value of each index; PC1 displays the first component of SPCA, which can aggregate the most information of five indicators; ECEI_{season_o} means the origin ECEI_season value; ECEI_{season_o_i}, ECEI_{season_o_min}, and ECEI_{season_o_max} represent the i_th, min and max value of ECEI_{season_o}; ECEI_season is the normalized value of ECEI_{season_o}; ECEI is the average value of all seasons in each period. ECEI ranges from 0 to 1, and a higher value indicates a higher EEQ.

2.3.2. Calculation of EI

Based on Technical Criterion for Ecosystem Status Evaluation (HJ 192-2015) [58], the calculation of EI consists of five parts, which are the biological abundance index (BAI), vegetation cover index (VCI), water network density index (WNDI), land stress index (LSI), and pollution load index (PLI). Equations (15)–(19) were used to calculate each index.

$$BAI=\frac{BI+HQI}{2}$$

(15)

where BI represents the biodiversity index, calculated from species diversity level, ecosystem diversity level, and landscape diversity level [59,60].

Table 2. Weight and sub-indicator system of BI.

Level	Weight	Sub-Indicator	Weight
Species diversity	0.60	Habitat quality index (HQI)	0.30
		Enhanced vegetation index (EVI)	0.15
		Water network density index (WNDI)	0.15
Ecosystem diversity	0.15	Percentage of habitat area (Sp)	0.05
		Simpson diversity index (SIDI)	0.10
Landscape diversity	0.25	Splitting index (SPLIT)	0.15
		Contagion index (CONTAG)	0.10

displays the weight and sub-indicator names for BI calculation. Among them, the Simpson diversity index (SIDI), Splitting index (SPLIT), and Contagion index (CONTAG) were calculated using Fragstats 4.2 software. HQI indicates habitat quality index, which is acquired by using InVEST model [61,62]. The weights of all indicators were determined by analytic hierarchy process method (AHP) [60].

$$NDV{I}_{avg}=A_{veg}\times \left(\frac{{\displaystyle {\sum }_{i=1}^n{P}_i}}{n}\right)$$

(16)

where P_i is the average value of the maximum monthly value of NDVI from May to September; n is the pixel number; A_veg is the normalization coefficient value; NDVI_avg is the value of VCI.

$$\mathit{WNDI}=\left[\begin{array}{c}{\mathit{A}}_{\mathit{riv}}\times {\mathit{Length}}_{river}{/\mathit{Area}}_{region}{+\mathit{A}}_{\mathit{lak}}\times Are{a}_{water}{/\mathit{Area}}_{region}\\ {+\mathit{A}}_{\mathit{res}}\times {\mathit{Souce}}_{\mathit{water}}{/\mathit{Area}}_{region}\end{array}\right]/3$$

(17)

where A_riv, A_lak, and A_res are the normalized coefficient values; Length_river is the length of all rivers in Zhaotong; Area_water is the area of water in Zhaotong, including lakes, reservoirs, rivers, canals, and offshore; Area_region is Zhaotong’s area; Source_water is the water resource amount [63].

$$LSI=A_{ero}\times \left(\begin{array}{l}0.4\times Are{a}_{severe}+0.2\times Are{a}_{middle}\\ +0.2\times Are{a}_{other}+0.2\times Are{a}_{building}\end{array}\right)/Are{a}_{region}$$

(18)

where A_ero is the normalized coefficient value; Area_severe, Area_middle, Area_other, and Area_building are the area of severe soil erosion, medium soil erosion, other soil erosion, and building land; Area_region is the city area.

$$\begin{array}{cc}PLI& =0.20\times {\mathit{A}}_{COD}\times {\mathit{Ems}}_{COD}/{\mathit{Pre}}_{annual}+0.20\times {\mathit{A}}_{N{H}_3}\times {\mathit{Ems}}_{N{H}_3}/{\mathit{Pre}}_{annual}+\\ & 0.20\times {\mathit{A}}_{S{O}_2}\times {\mathit{Ems}}_{S{O}_2}/Are{a}_{region}+0.10\times {\mathit{A}}_{YFC}\times {\mathit{Ems}}_{YFC}/Are{a}_{region}+\\ & 0.20\times {\mathit{A}}_{NOX}\times {\mathit{Ems}}_{NOX}/Are{a}_{region}+0.10\times {\mathit{A}}_{SOL}\times \mathit{V}{\mathit{ol}}_{SOL}/Are{a}_{region}\end{array}$$

(19)

where A_COD, A_NH3, A_SO2, A_YFC, A_NOX, and A_SOL are normalized coefficient values; Ems_COD, Ems_NH3, Ems_SO2, Ems_YFC, Ems_NOX, and Vol_SOL represent emissions of COD, NH₃, fume, and dust, NO_X and solid waste disposal volume; Pre_annual is the average annual precipitation; Area_region is the city area.

Due to the lack of long-term NH₃ and NO_X data, this study used CHAP datasets to represent the emission amount of NH₃ and NO_X [64,65]. As CHAP datasets represented spatial data, each county’s emission amount was allocated by the AHP method.

Table 3. Pollutant allocation indicator system.

Target Level	Element Layer	Indicator	Weight
Pollutant emissions allocation	Environment	Water quality (WQ)	0.38
	Economy	Gross domestic product (GRP)	0.33
	Society	Population density (PD)	0.29

shows the pollutant allocation indicator system.

2.3.3. Calculation of SDE

The SDE, long served as a general geographic information system tool, is a statistical method to measure the distribution characteristics of spatial features [66]. It can be constructed by ArcGIS directional distribution tool and can describe the characteristics of gravity center trend, dispersion, and directional trend [66]. In this study, SDE was used to present the aggregation characteristics of ECEI spatial distribution. One standard deviation was selected to perform SDE [67]. Equations (20)–(24) were used to calculate the gravity center, corner, and major and minor axes of SDE.

$$SD{E}_x=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(x_i-\overline{X}\right)}^2}}{n}}$$

(20)

$$SD{E}_y=\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left(y_i-\overline{Y}\right)}^2}}{n}}$$

(21)

$$\tan \theta =\frac{\left({\displaystyle \sum _{i=1}^n{\overline{x}}_i^2}-{\displaystyle \sum _{i=1}^n{\overline{y}}_i^2}\right)+\sqrt{{\left({\displaystyle \sum _{i=1}^n{\overline{x}}_i^2}-{\displaystyle \sum _{i=1}^n{\overline{y}}_i^2}\right)}^2+4{\left({\displaystyle \sum _{i=1}^n{\overline{x}}_i{\overline{y}}_i}\right)}^2}}{2{\displaystyle \sum _{i=1}^n{\overline{x}}_i{\overline{y}}_i}}$$

(22)

$${\sigma }_x=\sqrt{2}\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left({\overline{x}}_i\cos \theta -{\overline{y}}_i\sin \theta \right)}^2}}{n}}$$

(23)

$${\sigma }_y=\sqrt{2}\sqrt{\frac{{\displaystyle \sum _{i=1}^n{\left({\overline{x}}_i\sin \theta +{\overline{y}}_i\cos \theta \right)}^2}}{n}}$$

(24)

where x_i and y_i display the spatial location of each pixel; $\overline{X}$ and $\overline{Y}$ indicate the arithmetic mean center; ${\overline{x}}_i$ and ${\overline{y}}_i$ represent the difference between mean center and coordinate x and coordinate y; θ presents the corner. In this study, it represents the primary spatial distribution trend of ECEI.

2.3.4. Trend Analysis Method

Trend analysis is a method to analyze the change trend and intensity of a variable by performing a univariate linear regression analysis. The slope of the regression equation indicates the variable’s change trend [68,69]. In this study, trend analysis was applied to analyze the change trend of ECEI. Equation (25) shows the calculation formula.

$${\theta }_{slope}=\frac{n\times {\displaystyle \sum _{i=1}^{n}i\times N_i}-{\displaystyle \sum _{i=1}^{n}i}{\displaystyle \sum _{i=1}^n{N}_i}}{n\times {\displaystyle \sum _{i=1}^n{i}^2}-{\left({\displaystyle \sum _{i=1}^{n}i}\right)}^2}$$

(25)

where θ_slope denotes the changing trend; n is the number of study periods. In this study, n equals 5; N_i is the ith value of ECEI; θ_slope > 0 indicates an increasing trend, and vice versa.

3. Results

3.1. Spatiotemporal Change Analysis of ECEI

Figure 3 shows the average ECEI values of Zhaotong from 2000 to 2020. Based on Figure 3, all ECEI values in Zhaotong were higher than 0.690, indicating an inverted “N” change trend. In 2005, the ECEI had the lowest value (0.693). In 2015, the ECEI had the highest value (0.749). Generally, the EEQ of Zhaotong has shown an upward trend over the past 20 years, with the ECEI value increasing by 0.23%.

To understand the spatial characteristics of the ECEI, a spatial distribution map at grid and county levels was drawn (Figure 4 and Figure 5).

According to Figure 4, the distribution pattern of EEQ in Zhaotong was relatively stable. Regions located in the northern and eastern parts had high ECEI values. Western areas had moderate ECEI values, better contiguity, and rich patch detail. Areas with low ECEI values showed highly fragmented characteristics. In addition, these regions were mainly distributed in Weixin county, Zhenxiong county, Zhaoyang district, surrounding areas of Shuifu city, Suijiang county, Yongshan county, Daguan county, and Qiaojia county. Except for Chaoyang district, all remaining regions displayed an upward trend to varying degrees. The EEQ of the southern Zhaoyang district was consistent in the past 20 years, with no expansion. In 2015, areas with low ECEI values in Zhaoyang district showed a slight reduction. In addition, the EEQ in the western regions, including Yongshan county, Zhaoyang district, Ludian county, and Qiaojia county, have improved. Due to the lack of satisfying images, in 2015, some regions did not have valid ECEI values. Here, to fully understand the area change with respect to different grades, we set 0.2 as the interval to equally divide ECEI into five grades, which were bad [0,0.2], fair (0.2,0.4], moderate (0.4,0.6], good (0.6,0.8] and excellent (0.8,1].

Table 4. Area percentage of different ECEI grades.

ECEI Grades	2000	2005	2010	2015	2020
	Area Percentage/%	Area Percentage/%	Area Percentage/%	Area Percentage/%	Area Percentage/%
Bad	0.37	0.49	0.61	0.79	0.73
Fair	13.45	13.48	13.48	13.86	10.04
Moderate	22.39	22.31	22.36	22.44	23.70
Good	0.57	2.65	0.59	0.63	7.27
Excellent	63.22	61.07	62.96	62.28	58.26

shows the area percentage of different ECEI grades. Figure 6 shows the area percentage bar chart of Zhaotong from 2000 to 2020.

Data from Table 4 and Figure 6 show that the area percentage of different grades in different years was relatively stable. Specifically, for the bad grade, its area percentage displayed an upward trend, with its value increasing from 0.37% in 2000 to 0.73% in 2020. However, the percentage of this grade in all periods was small (lower than 1%). For the fair grade, its area percentage presented a downward change trend, especially in 2020; the area percentage decreased from 13.86% in 2015 to 10.04% in 2020. For the moderate grade, its area percentage increased from 22.39% in 2000 to 23.70% in 2020. For the good grade, its area percentage showed a fluctuating upward trend, with the area percentage rising from 0.57% in 2000 to 7.27% in 2020. As for the excellent grade, its area percentage displayed a fluctuating downward change trend, with the area percentage decreasing from 63.22% in 2000 to 58.26% in 2020. Generally, in 2020, the area percentage of excellent grade decreased. However, the area percentage of moderate, good, and excellent grades in 2020 was the highest. In addition, in all periods, the area percentage of the excellent grade was the highest, indicating that the overall EEQ of Zhaotong city was good.

3.2. Change Direction of ECEI Based on SDE

Figure 7 shows the change direction of ECEI in Zhaotong from 2000 to 2020.

According to Figure 7, the θ range of SDE was from 57.06° to 62.90° with the central areas almost located in the SDE. The distribution pattern of EEQ in Zhaotong presented the direction of “southwest–northeast.” Based on the standard deviation change situation of the major and minor axes of the ellipse, this study was divided into four periods, which were 2000–2005, 2005–2010, 2010–2015, and 2015–2020. Specifically, from 2000 to 2005, the standard deviation of the major and minor axes of the ellipse increased, and regions with high ECEI values moved along the major axis and toward the northeast of Zhaotong city. At the same time, these regions spread toward the minor axis. From 2005 to 2010, the standard deviation of the major and minor axes of the ellipse indicated that the spatial aggregation situation of the ECEI moved toward the gravity center of the ellipse with the coordinate point changing from (386983.50, 3045862.44) to (381051.81, 3041515.11). In general, compared with the previous period, in this period, the changes in southwestern regions’ EEQ played an important role. From 2010 to 2015, standard deviation increments of the major and minor axes of the ellipse were small, and the SDE spatial distribution over two years was similar, indicating that the spatial distribution of EEQ was stable. From 2015 to 2020, standard deviation increments of the major and minor axes of the ellipse were accelerated. Compared with the first period, the transport and expansion intensity were stronger, indicating that the changes in the northeastern regions’ EEQ played an enhanced driving effect. In addition, in this period, the gravity center was located in the northeast corner of Zhaoyang district. In general, from 2000 to 2020, the standard deviation of the major axis of the ellipse moved toward the northeast of Zhaotong city with θ of SDE changing from 57.06° to 62.90°, indicating the improvement of northeastern regions’ EEQ.

3.3. Change Trend of ECEI Based on the Trend Analysis Method

Figure 8 shows the spatial distribution of the ECEI changing trend in Zhaotong from 2000 to 2020.

According to Figure 8, the change trend value ranged from −0.24 to 0.24. Change trend values of most regions were close to 0, indicating that the EEQ of those regions did not show a noticeable improvement or decrease. Regions with low change trend values were mainly distributed in the center of Zhaoyang district. These regions primarily comprised artificial surfaces with intense human activities. To further understand these changes, the center regions of Zhaoyang district were enlarged (b). In Figure 8b, regions with the red color displayed an obvious eco-environmental deterioration. Combined with Figure 8c,d, these regions belonged to urban-expanded areas. In 2021, the urbanization rate of Zhaotong city reached 40.50%, and the build-up area of Zhaoyang district was 62 km². Other areas with low change trend values presented a scattered distribution. Regions with high change trend values also displayed a scattered distribution. Generally, in the past 20 years, with the increase in urbanization levels, especially the increase in artificial surfaces, the EEQ of the center of Zhaoyang district displayed an evident decrease trend. However, the EEQ’s change trend in most regions of Zhaotong city remained stable, with the change trend value close to 0. Other regions with high or low change trend values showed a scattered distribution.

4. Discussion

4.1. Analysis of ECEI’s SPCA Results

Based on the aforementioned analysis, each ECEI season’s SPCA results were analyzed. Among them, the eigenvalue and percentage of PC1 are shown in

Table 5. Each ECEI season’s SPCA results.

Year	Season	Eigenvalue of PC1	Percentage of PC1/%
2000	Spring–Summer	0.0303	73.40
	Autumn–Winter	0.0292	80.24
2005	Spring–Summer	0.0325	69.94
	Autumn–Winter	0.0298	73.72
2010	Spring	0.0312	72.16
	Summer	0.0302	70.92
	Autumn	0.0303	65.83
	Winter	0.0325	69.94
2015	Spring	0.0321	81.47
	Summer	0.0290	86.93
	Autumn	0.0297	85.43
	Winter	0.0308	78.53
2020	Spring	0.0352	77.91
	Summer	0.0306	80.04
	Autumn	0.0301	82.77
	Winter	0.0313	79.18

.

According to Table 5, in 2000 and 2005, due to the lack of cloudy-free images, only two seasons were taken into consideration. The eigenvalue represented the information aggregated by the SPCA method. A higher eigenvalue aggregated more information. However, there was no comparison between two or more periods. For example, for the spring–summer season in 2000, the eigenvalue of PC1 was the highest compared with PC2 to PC5. Specifically, for the three remaining periods, four seasons’ cloudy-free images were all available. Specifically, in 2000, the two seasons’ eigenvalues of PC1 were 0.0303 and 0.0292, with PC1 percentages of 73.40% and 80.24%. In 2005, the two seasons’ eigenvalues of PC1 were 0.0325 and 0.0298, with the PC1 percentages of 69.94% and 73.72%. In 2010, four seasons’ eigenvalues of PC1 were 0.0312, 0.0302, 0.0303, and 0.0325, with PC1 percentages of 72.16%, 70.92, 65.83%, and 69.94%. In 2015, four seasons’ eigenvalues of PC1 were 0.0321, 0.0290, 0.0297, and 0.0308, with percentages of 81.47%, 86.93%, 85.43%, and 78.53%. In the year 2020, the eigenvalues of different seasons were 0.0352, 0.0306, 0.0301, and 0.0313, and percentages of PC1 in four seasons were 77.91%, 80.04%, 82.77%, and 79.18%. Generally, the percentages of PC1 were all higher than 65%, with an average percentage value of 76.77%. These results were consistent with the existing RSEI studies. In these studies, even the results of the contribution of PC1 varied. Typically, the contribution was higher than 60% [70,71]. Therefore, based on Table 5, the first component of the SPCA could aggregate the most information of the five indicators, demonstrating the first component’s ability to influence the regional EEQ’s assessment.

4.2. Validation of ECEI by Comparing with EI

To further validate the effectiveness of the ECEI, the EI was adopted to compare with the ECEI. According to the aforementioned calculation steps, the results of the EI in 2010, 2015, and 2020 were calculated considering the data accessibility. Figure 9 shows the eight sub-indexes calculation results of BI in 2020.

Based on Figure 9, different sub-indexes showed noticeable differences. For WNDI, southern and northwestern Zhaotong had high WNDI values while central regions had low WNDI values, indicating that those regions’ water networks were relatively low. Areas with low HQI values were mainly distributed in Zhaoyang district as those regions had intense human activities. However, the HQI value in northern Zhaotong was relatively high. This is mainly because of those regions’ high vegetation coverage. The spatial distribution of the EVI further validated that areas with high vegetation index are mainly distributed in northern Zhaotong. For SPLIT, regions with high SPLIT values were primarily distributed in southwestern and western parts, indicating that landscape segmentation degrees of these regions were high. However, at the center of Zhaoyang district, the landscape segmentation degree was low. Combined with CONTAG and SIDI, those regions had contiguous artificial surfaces, which further proved that regions with intense human activities usually had low ecosystem and landscape diversity [72,73].

Figure 10 shows the spatial distribution map with the average EI value of Zhaotong city during three periods.

Based on Figure 10, the spatial distribution situation of the average EI during the three periods was similar. Specifically, the average EI values of northern Zhaotong were high, especially in Yanjin and Suijiang counties. However, Zhaoyang district and Ludian county had low average EI values. From the perspective of spatiality, average EI values presented the pattern of high values in the north and low values in the south. From the standpoint of the highest–lowest value, from 2010 to 2020, the lowest values were 60.08, 60.85, and 60.49, respectively, while the highest values were 71.85, 72.15, and 74.57. Generally, the lowest value showed a fluctuating upward trend, with the value increasing by 0.69%. The highest value showed an increasing upward trend with the value increased by 3.78%. These results indicated that the EEQ of Zhaotong improved from 2010 to 2020, which was consistent with the results of the ECEI. In addition, based on the technique criteria, all years’ highest and lowest EI values were higher than 55, indicating that all regions of Zhaotong city qualified for good grades.

To further validate the accuracy of the ECEI, average ECEI and RSEI values of the same three periods were calculated to draw the scatter diagram and establish the fitting curve (Figure 11).

Based on Figure 11, the coefficients of determination (R²) between the EI and ECEI were 0.95 (2010), 0.88 (2015), and 0.89 (2020), with an average R² of 0.90. As for the R² between the EI and RSEI, the values of the three periods were 0.31 (2010), 0.92 (2015), and 0.82 (2020), with an average value of 0.69. It could be found that the average R² value between the EI and ECEI was higher than that of the EI and RSEI, indicating that the ECEI could provide more accurate and scientific EEQ results when considering land surface and season differences. In addition, the R² of the ECEI was more stable than the RSEI with the value ranging from 0.88 to 0.95, while the R² of the RSEI went from 0.31 to 0.92. In addition, the R² of considering all values of the three periods further validated the accuracy of the ECEI with the value of 0.84, which was higher than that of the RSEI (R² = 0.56). Generally, after considering land surface and season differences, the evaluation of the EEQ in typical cloudy regions was more accurate and scientific compared with the RSEI.

4.3. Implications, Limitations, and Further Study

To evaluate the EEQ of cloudy regions, land surface and season differences were considered to construct a new index-ECEI to assess Zhaotong’s EEQ from 2000 to 2020 accurately and scientifically. Then, the ECEI and RSEI results were compared to validate the effectiveness of the ECEI. Based on the aforementioned analysis, Zhaotong’s EEQ showed a fluctuating upward trend in the past 20 years, with the value increasing from 0.729 to 0.730. From 2000 to 2005, ECEI’s value decreased from 0.729 to 0.693. In this period, due to the rapid population growth and low economic development level, people decided to get more food by destroying forests, cutting trees, and abandoning grazing, thus causing serious water and soil loss and forest coverage loss [72]. From 2002 to 2004, the area of water and soil loss in Zhaotong city increased from 11307.93 km² to 13362.49 km² [72,73]. From 2005 to 2020, with the implementation of a series policies, including the “Returning Farmland to Forest (grass) Project”, the forest coverage rate increased from 32.60% in 2011 to 47.80% in 2021. As a result, northern counties had a higher forest coverage rate compared with southern counties. Specifically, Suijiang county’s forest coverage rate reached 71.99% in 2021. In addition, Yanjin, Daguan, Yongshan, and Suijiang counties were listed in the National Key Ecological Function Areas, implemented by the Chinese government since 2007. The number of nature reserves at all levels increased from 8 in 2011 to 23 in 2021, with the area increasing to 1656.00 km². The controlled area of rocky desertification reached 1180 km². Moreover, the wastewater treatment rate and garbage harmless treatment rate increased from 41.80% and 45.80% to 98.40% and 100.00%, respectively. Generally, the spatial distribution of the EEQ displayed a pattern of high values in the north and low values in the south, which was consistent with the vegetation coverage situation. Based on the statistical yearbook of Zhaotong, at the end of 2021, the forest coverage rate of Suijiang county was the highest, with a value of 71.99%. Except for Suijiang county and other northern counties in Zhaotong city, forest coverage rates of Yanjin county and Shuifu county were 64.62% and 66.5%. However, the Zhaoyang district’s forest coverage rate was only 38.67%. These results indicated the accuracy of ECEI evaluation results. In addition, the EEQ value of Zhaoyang district was the lowest due to its intense human activities. Existing studies found that human activities could exert a significant influence on the regional EEQ [74,75]. In the future, regions with intense human activities, such as Zhaoyang district, should seek coordinated development between the economy and eco-environment [23].

This paper established the ECEI and evaluated the EEQ of Zhaotong in the past 20 years. Compared with the RSEI, the ECEI results were more accurate and stable. However, this paper still had some limitations. On the one hand, the regional eco-environment belongs to a compound system, including society, economy, and nature [23]. The RSEI only considers four main aspects of the eco-environment. In this study, land surface and season differences were taken into consideration. However, some components, such as air and soil environments, should also be taken into account in the evaluation model. Hence, in further studies, relevant air and soil indicators should be carefully selected to establish a more comprehensive index. On the other hand, for cloudy regions, although the adjacent year was also used to select the satisfying images, this study found that some places did not have any satisfying images. Hence, in further studies, image fusion and image construction from multi-satellites may provide helpful ways to evaluate the regional EEQ thoroughly. Generally, despite these shortcomings, this study still provided a new framework for assessing cloudy regions’ EEQ more accurately and scientifically.

5. Conclusions

Based on multi-source datasets, this study constructed a new eco-environmental comprehensive evaluation index (ECEI) by considering land surface and season differences to evaluate Zhaotong’s EEQ in the past 20 years. Five aspects, such as the biological abundance index (BAI), vegetation cover index (VCI), water network density index (WNDI), land stress index (LSI), and pollution load index (PLI), were also calculated to acquire the ecological index (EI) of Zhaotong at the county level from 2010 to 2020. Finally, spatiotemporal characteristics, standard deviation ellipse (SDE), and trend analysis of Zhaotong’s EEQ at different periods were analyzed. In addition, scatter diagrams of EI–ECEI and EI–RSEI were explored. The main conclusions were as follows: (1) The EEQ of Zhaotong showed a fluctuating upwards trend (0.0058 a⁻¹), with average ECEI values of 0.729, 0.693, 0.722, 0.749, and 0.730. Since 2005, a series of policies, including the “Returning Farmland to Forest (grass) Project”, were implemented to improve Zhaotong’s EEQ. (2) The spatial distribution pattern of EEQ showed high values in the north and low values in the south, which was consistent with the actual forest coverage rate distribution. For example, Suijiang county had a forest coverage rate of 71.99% in 2021. Zhaoyang district had the lowest ECEI value, which was mainly caused by urban sprawl and intense human activities. (3) From 2000 to 2020, the standard deviation of the major axis of the ellipse moved northeast of Zhaotong city with θ of SDE changing from 57.06° to 62.90°, indicating the improvement of northeastern regions’ EEQ; this was related to the increase of National Key Ecological Function Area and nature reserves numbers in northern counties. (4) The coefficients of the determinant (R²) between the EI and ECEI were 0.84, which was higher than that of EI–RSEI (R² = 0.56). Generally, the ECEI could provide more accurate spatial EEQ evaluation results compared with the RSEI, which could contribute to the regional eco-environmental management, as the spatial pattern of the regional EEQ could be quickly obtained and actual measures could improve the regional EEQ with low ECEI values.