Urban ecosystem quality assessment based on the improved remote sensing ecological index

Study area

The municipal area in Chongqing covers an area of 82,400 km² and is divided into three main functional zones: the city proper, the Three Gorges Reservoir area in the northeast, and the Wuling Mountain area in the southeast. Chongqing’s city proper comprises 21 districts, including the central urban area (Yuzhong, Dadukou, Jiangbei, Nan’an, Shapingba, Jiulongpo, Beibei, Yubei, and Ba’nan) and the main urban area (Fuling, Changshou, Jiangjin, Hechuan, Yongchuan, Nanchuan, Qijiang, Dazu, Bishan, Tongliang, Tongnan, and Rongchang), covering a total area of 28,700 km² (Fig. 1). Accounting for 35% of Chongqing’s land area, this region accommodates 66% of the city’s permanent residents and 74% of the city’s urban population, and it contributes 77% of the regional GDP. It is characterized by a relatively high level of urbanization and is closely linked to urban development.

Data and preprocessing

In this study, data for the city proper in Chongqing were collected for the years 2002, 2006, 2010, 2014, 2018, and 2022. Table 1 summarizes the key data used in this study (land cover data and statistical data for calculating the EI values). Due to the similarity of the annual trends of the RSEI and the average summer RSEI values, the data for June to August in each target year were selected (Ji et al., 2022; Zhang et al., 2021).

All image processing was performed on Google Earth Engine (GEE) using Landsat 5 (LC05) and Landsat 8 (LC08) data. The specific dataset identifiers were: LANDSAT/LT05/C02/T1_L2 for Landsat 5 (EROS, 2020a), and LANDSAT/LC08/C02/T1_L2 for Landsat 8 (EROS, 2020b). Prior to release, both satellite datasets underwent radiometric calibration and atmospheric correction. These preprocessing steps effectively eliminated sensor response variations, temporal imaging differences, and atmospheric interference, ensuring accurate representation of surface reflectance. LC05 and LC08 images were calibrated from sensor radiance to surface reflectance using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al., 2006) and the Landsat Surface Reflectance Code (LaSRC) (Vörösmarty, Pahl-Wostl & Bhaduri, 2013), respectively. These datasets underwent preprocessing steps.

Due to Chongqing’s humid and rainy climate, the images were often affected by clouds. To address cloud contamination, we implemented the Fmask algorithm for cloud removal (Foga et al., 2017; Mateo-García et al., 2018; Zhu & Woodcock, 2012). In GEE, Fmask is encapsulated as a convenient function, allowing cloud and cloud shadow masking by extracting the “StateQA” band. After cloud removal, the images were composited using the median value. The effectiveness of cloud removal was evaluated through cloud cover calculations and visual analysis. Specifically, the “QA_PIXEL” band was used to identify cloudy pixels, and cloud coverage was calculated as the ratio of cloudy pixels to the total number of pixels. The results showed that after cloud removal, the cloud coverage of the six composite images from 2002 to 2022 was close to zero, meeting the requirements for further analysis. Since the cloud removal method strictly filters out cloud-contaminated areas, regions that are severely affected by cloud contamination in all available images may lose valid data after cloud removal, resulting in missing areas. Therefore, after cloud removal, it is necessary to fill the missing areas using images from the same month of the previous or following year.

Calculation of various indicators

(1) NDVI: The NDVI has been widely used to assess and monitor the health and coverage of surface vegetation. It is calculated as follows: (1)${\displaystyle NDVI=\frac{B_{Nir}-B_{Red}}{B_{Nir}+B_{Red}}}$

where B_Nir is the reflectance in the near–infrared band, and B_Red is the reflectance in the red band.

(2) WET: The moisture index not only represents open water bodies but is also closely related to the moisture in the soil and vegetation. Before extracting the WET index, water bodies are masked using the modified normalized difference water index (Xiong et al., 2021) (MNDWI, Eq. (3)). The formula for calculating the WET index (Crist, 1985) is ${\displaystyle WE{T}_{TM}=0.0315B_{Blue}+0.2021B_{Green}+0.3102B_{Red}}$ ${\displaystyle +0.1594B_{Nir}-0.6806B_{Swir1}-0.6109B_{Swir2}}$ ${\displaystyle WE{T}_{OLI}=0.1511B_{Blue}+0.1973B_{Green}+0.3283B_{Red}}$ (2)${\displaystyle +0.3407B_{Nir}-0.7117B_{Swir1}-0.4559B_{Swir2}}$ (3)${\displaystyle MNDWI=\frac{B_{Green}-B_{Swir1}}{B_{Green}+B_{Swir1}}}$

where B_Blue, B_Green, B_Red, B_Nir,  B_Swir1, and B_Swir2 are the reflectances of the Landsat data in the blue, green, red, near–infrared, shortwave infrared 1, and shortwave infrared 2 bands, respectively.

(3) LST: The LST represents the surface temperature in an area. In this study, the Landsat LST product (Ermida et al., 2020) provided by the GEE platform was utilized, requiring only conversion to Celsius. The calculation formula is as follows: (4)${\displaystyle LST=B_{ST}-273.15}$

where B_ST is the land surface temperature band (for the LANDSAT/LT05/C02/T1_L2 dataset, it is the ST_B6 band, and for the LANDSAT/LC08/C02/T1_L2 dataset, it is the ST_B10 band). To convert the temperature to Celsius, 273.15 is subtracted from the Kelvin value.

(4) NDSI: The NDSI consists of two components: the scarcity index (SI), which reflects the extent of vegetated land scarcity (Rikimaru, Roy & Miyatake, 2002); and the index-based built-up index (IBI), which reflects built-up land conditions (Xu, 2008). Its calculation formula is

${\displaystyle NDSI=\frac{SI+IBI}{2}}$ ${\displaystyle SI=\frac{\left(B_{Swir1}+B_{Red}\right)-\left(B_{Blue}+B_{Nir}\right)}{\left(B_{Swir1}+B_{Red}\right)+\left(B_{Blue}+B_{Nir}\right)}}$ (5)${\displaystyle IBI=\frac{\frac{2B_{Swir1}}{B_{Swir1}+B_{Nir}}-\left[\frac{B_{Nir}}{B_{Red}+B_{Nir}}+\frac{B_{Green}}{B_{Swir1}+B_{Green}}\right]}{\frac{2B_{Swir1}}{B_{Swir1}+B_{Nir}}+\left[\frac{B_{Nir}}{B_{Red}+B_{Nir}}+\frac{B_{Green}}{B_{Swir1}+B_{Green}}\right]}}$

where B_Blue, B_Green, B_Red, B_Nir, B_Swir1, and B_Swir2 are the reflectances of the Landsat data in the blue, green, red, near–infrared, shortwave infrared 1, and shortwave infrared 2 bands, respectively.

(5) DI: Feng, Feng & Feng (2018) developed the particle difference index (DI) to represent changes in particle concentration based on the characteristics of the PM_2.5, which increases the reflectance in the red band and decreases the reflectance in the near–infrared band. Its calculation formula is (6)${\displaystyle DI=B_{Red}-B_{Nir}}$

where B_Red and B_Nir are the apparent reflectance or radiance in the red and near–infrared bands, respectively.

DRSEI comprehensive index development

The four indicators included in the RSEI are greenness, moisture, dryness, and heat. Greenness is represented by the normalized difference vegetation index (NDVI), moisture is represented by the humidity index obtained through remote sensing tasseled cap transformation (WET), dryness is represented by the average of the soil index (SI) and the built–up index (IBI) (NDSI), and heat is represented by the land surface temperature (LST) derived from Landsat. Thus, the definition of the RSEI is (Xu, 2013a). ${\displaystyle RSEI=f\left(NDVI,WET,LST,NDSI\right).}$

By replacing LST in the RSEI with the difference index (DI), a new type of remote sensing ecological index is constructed, denoted as the DRSEI. Therefore, the definition of the DRSEI is ${\displaystyle DRSEI=f\left(NDVI,WET,NDSI,DI\right).}$

To integrate these four indicators into a single index, the principal component analysis (PCA) method (Bro & Smilde, 2014) is commonly employed. PCA identifies the correlations among these indicators by analyzing their variations and determines their respective weights based on each indicator’s contribution to the principal components. This approach transforms the original data into a few representative principal components, thereby reducing data complexity and better capturing the overall ecological environment. Therefore, the formula for the DRSEI is: (7)${\displaystyle DRSEI=PC1\left[f\left(NDVI,WET,NDSI,DI\right)\right]}$

where PC1 represents the first principal component obtained from PCA.

For the calculated four indicators, due to their non–uniform dimensions, it was necessary to normalize these indicators before PCA (Ye & Kuang, 2022). That is, their values needed to be normalized to the range of [0,1]. The normalization formula is (8)${\displaystyle N_i=\frac{I_i-I_{min}}{I_{max}-I_{min}}}$

where N_i is the normalized value of a certain indicator, I_i is the value of that indicator in pixel i, and I_min and I_max are the minimum and maximum values of that indicator, respectively.

PCA was performed on the normalized indicators to retain the primary information in the data while achieving reduction of the dimensionality. To enable measurement and classification of the DRSEI, it was also normalized. The formula for this normalization is (9)${\displaystyle DRSEI=\frac{DRSE{I}_i-DRSE{I}_{min}}{DRSE{I}_{max}-DRSE{I}_{min}}}$

where DRSEI_i is the value of DRSEI for pixel i, and DRSEI_min and DRSEI_max are the minimum and maximum values of DRSEI, respectively.

Model verification

The optimization effect of DRSEI, as an improved version of RSEI, can be validated from two dimensions:

1. Accuracy improvement assessment: using the ecological index (EI) as the evaluation standard, the improvement of DRSEI in ecological quality assessment is quantitatively evaluated by comparing the degree of convergence between DRSEI, RSEI, and EI. The selection of EI as the reference benchmark is based on its widespread adoption in China’s environmental statistical yearbooks for ecological quality assessment. The fit between DRSEI, RSEI, and EI is examined using three core indicators, with specific calculation formulas provided in Table 2: correlation coefficient, root–mean–square difference, and standard deviation.

2. Regional applicability analysis: a multi-regional comparative experiment is conducted to analyze the applicability differences between DRSEI and RSEI in regions with distinct ecological characteristics, assessing their stability and generalizability in complex geographical environments. Specifically, the DRSEI and RSEI scores for different ecological regions are calculated and systematically compared with the EI reference values for these regions to reveal the adaptability of the models under varying geographical conditions.

Comparison of DRSEI and RSEI in the city proper in Chongqing

PCA was implemented on the GEE platform (Google Earth Engine) to obtain the loading coefficients of each indicator and the contribution of the first principal component (Table 3). The results indicate that the first principal component of the DRSEI accounted for over 75% of the total variance, demonstrating a capacity equivalent to that of the RSEI in synthesizing multi-indicator information. Compared to the RSEI, the higher variance contribution of the first principal component in the DRSEI suggested that its model accounted for a larger proportion of the total variance. This implied that the DRSEI summarized the data’s main characteristics more effectively than the RSEI. Consequently, it enhanced the capability to characterize ecological factor trends. Regarding loading polarities, greenness and wetness showed positive loadings, indicating beneficial contributions to the ecosystem, while dryness and air pollution indices exhibited negative loadings, reflecting detrimental impacts on EQ. These findings aligned closely with the RSEI pattern. Additionally, across all study years, the consistent signs of indicator loadings in the DRSEI suggested that the model had a stable capacity to reflect EQ evolution.

Both the contribution rate of the first principal component and PCA-derived indicator loadings suggested that the DRSEI enhanced multi-indicator integration; however, its validity for EQ evaluation required further verification. We selected 2018 and 2022 as representative cases where both indices showed high contribution rates, then conducted comparative analyses between the indices and their corresponding EI values across administrative divisions.

At the county level, we computed correlation coefficients, RMS differences, and standard deviations between EI and both RSEI/DRSEI for 2018 and 2022. The correlation coefficient revealed which model better captured the authentic ecological variation trends overall. RMS differences quantified deviations between EI and the indices (EI-RSEI/DSEI pairs), while standard deviations measured the magnitude of data dispersion, thereby reflecting the spatial heterogeneity of ecological indices. To more intuitively illustrate the relationships among the three evaluation indices, we used a Taylor diagram (Fig. 2). The Taylor diagram can simultaneously display the correlation, standard deviation, and RMS difference of multiple models relative to the reference value (EI) within the same coordinate system, overcoming the limitations of traditional tables or single scatter plots in comprehensively comparing multiple indices. Each point in the diagram represents a model. The angle between a radial axis and the abscissa represents the correlation coefficient. Smaller angles (approaching 0°) indicate values closer to 1, denoting stronger linear consistency with ground truth data. The distances along orthogonal axes denote standard deviations. Euclidean distances between model markers and the reference point represent RMS differences; shorter distances indicate smaller RMS errors (Taylor, 2001).

First, we compared the standard deviations in the Taylor diagram. Analysis of the relative positions between scatter points and the solid arc connected to the REF point in Fig. 2 revealed that standard deviations followed the order: RSEI >DRSEI >EI. However, the differences among them were relatively small, indicating that while both DRSEI and RSEI exhibited slightly greater data variability than EI, the extent of this difference was limited. This suggests that these models may slightly amplify the magnitude of ecological changes. Nevertheless, the comparison demonstrates the DRSEI’s standard deviation is closer to EI’s value, suggesting its superior realism in EQ representation. Next, we compared the correlation coefficients, which correspond to the angles of the radial lines in the diagram. The results showed that both indices had correlation coefficients of approximately 0.6 with the EI, indicating moderate correlations. Notably, the DRSEI exhibited a higher correlation coefficient than the RSEI, which suggests that the spatial distribution of the DRSEI was more synchronized with the actual EQ represented by the EI. Finally, we examined the root-mean-square (RMS) difference, which quantifies the overall discrepancy between the model values and the actual measurements. This metric is visualized in the diagram as the distance between scatter points and the reference (REF) point. Comparison of scatter point positions relative to the pink arc shows the DRSEI exhibiting smaller RMS difference, demonstrating superior fitting accuracy and spatial consistency in EQ assessment. These findings further validated the improvements in the DRSEI.

Compared to the RSEI, the DRSEI exhibited two key improvements in EQ assessment: fluctuations closer to true values; and higher spatial synchronization coupled with smaller overall error. These findings demonstrated its improved accuracy.

Comparison of DRSEI and RSEI in other cities

The DRSEI performed well in the main urban area of Chongqing, but its effectiveness might vary across cities with distinct geographical locations and climatic conditions. To systematically compare the DRSEI and the RSEI, this study selected Beijing (a northern city) and Guangzhou (a southern city) for comparison. Beijing is located in the northern North China Plain, while Guangzhou lies in the coastal area of South China. These significant latitudinal and geographical contrasts provided an ideal foundation for assessing the applicability of the indices across diverse environments. For Beijing, district-level ecological index (EI) data for 2022 were obtained from the Beijing Ecological Environment Status Bulletin. In contrast, Guangzhou’s EI reporting ceased after 2020, requiring the use of 2020 data. Differences between DRSEI, RSEI and EI were calculated and visualized using lollipop plots, explicitly demonstrating their discrepancies, as shown in Fig. 3. The lollipop plot not only intuitively displays the magnitude of the deviations of DRSEI and RSEI from EI across districts but also clearly indicates the direction of these deviations.

Figure 3 reveals that both DRSEI and RSEI show limited effectiveness in evaluating Beijing’s district-level EQ. Both models demonstrate systematic underestimation across the municipality. However, three key observations emerge: First, discrepancies exist between citywide and district-level evaluations. For both models, citywide EQ assessments exhibited higher accurate compared to district-level evaluations. Second, at the district level, the DRSEI achieved higher accuracy in central urban areas (including Dongcheng, Xicheng, Chaoyang, Fengtai, Shijingshan, Haidian, Mentougou, and Daxing). Smaller discrepancies between DRSEI and EI values in these districts demonstrate its enhanced suitability for urbanization-dominated ecosystems. In contrast, the RSEI demonstrated higher accuracy in assessing Beijing’s suburban districts (Fangshan, Shunyi, Huairou, Pinggu, Miyun, and Yanqing). These areas typically feature higher vegetation cover and lower urbanization levels, suggesting that the RSEI may be more suitable for reflecting EQ in such natural environments.

However, Guangzhou’s EQ assessment revealed marked disparities between DRSEI and RSEI performances. The RSEI generally overestimated EQ, while the DRSEI consistently underestimated it. This systematic divergence may stem from differences in model sensitivity to ecological factors under the humid and vegetation-rich climatic of southern China. Compared to district-level evaluations, the DRSEI demonstrated greater accuracy in citywide EQ assessments. This again demonstrates that, at a large scale, the DRSEI model can more stably capture the overall trend of ecological change. Additionally, the DRSEI showed significantly higher accuracy than the RSEI in assessing Guangzhou’s central urban districts, including Liwan, Yuexiu, Haizhu, Tianhe, and Baiyun. This spatial consistency with Beijing’s pattern strengthens evidence for differential model applicability across urbanization gradients. The DRSEI is better suited for assessing EQ in highly urbanized areas, while the RSEI performs better in natural ecological environments.

As DRSEI and RSEI demonstrate distinct regional applicability advantages and diverge in key metrics (DI for DRSEI versus LST for RSEI), systematically investigating LST and DI as potential core drivers becomes critical for elucidating the mechanistic basis of model performance variations. To validate LST and DI impacts on EQ, we first analyzed correlations among these parameters and EI (Fig. 4). The correlation analysis aims to assess whether LST and DI exhibit equivalent significance for ecological quality, thereby investigating their roles in the performance discrepancies between models. Should DI demonstrate a statistically stronger correlation with EI compared to LST, this would imply that the superiority of DRSEI in ecological quality assessment may originate from the validity of DI as a surrogate indicator. Specifically, DI could provide more precise characterization of ecological conditions, whereas LST might play a comparatively subordinate role. Conversely, if LST and DI show comparable correlations with EI, it would suggest both parameters exert similar influences on ecological quality. In such cases, the accuracy differences between DRSEI and RSEI might not be exclusively attributed to the substitution of evaluation factors, but could also be attributable to additional confounding variables.

Beijing showed DI-EI and LST-EI correlation coefficients of −0.92 (R² = 0.86) and −0.93 (R² = 0.87), respectively. Guangzhou demonstrated stronger DI-EI (−0.96, R² = 0.93) and LST-EI (−0.96, R² = 0.92) correlations. The consistently strong negative correlations (high R² values) between DI/LST and EI across both cities confirm their significant influence on EQ. These robust correlations substantiate the rationale and feasibility of substituting the DI for the LST in EQ assessments. Notably, despite the strong correlations between the DI and the LST with the EI in both cities, the overall evaluation results of the DRSEI and the RSEI in Guangzhou showed substantial differences. This suggested that the discrepancy in evaluation accuracy was not entirely caused by the fundamental relationship between the LST and the DI with EQ.

The PCA weighting mechanism is influenced by the spatial variability of indicators: regions with greater variation exhibit higher variance contributions and loading weights. This necessitates examining how indicators’ spatial distribution characteristics interact with the model’s weighting mechanism (Zheng et al., 2022). From a geospatial perspective, this study aimed to analyze the formation path of model differences. To elucidate this mechanism comprehensively, the study integrates analysis of spatial frequency distributions and principal component loadings for both DRSEI and RSEI across Guangzhou and Beijing. The spatial frequency histograms (Fig. 5) characterize the value distribution patterns of individual indicators across distinct geographical zones, particularly revealing whether specific indicators demonstrate localized concentration or extensive dispersion. The principal component loadings (Table 4), conversely, quantify each indicator’s contribution during the PCA weighting procedure. By synthesizing these spatial distribution patterns with their corresponding statistical weights, this investigation systematically demonstrates how spatial heterogeneity of indicators modulates PCA weighting schemes, and further reveals how such modulations generate regional discrepancies in model applicability between DRSEI and RSEI. These findings collectively advance our understanding of the mechanistic differences underlying the two indices.

Since the NDVI, WET, and NDSI are included in both models, this study focused on the differences between the LST and the DI. As shown in Fig. 5, LST and DI exhibit significant spatial variations in Beijing; their data distribution patterns show notable similarity. When combined with the overall scores of the DRSEI and the RSEI in Beijing shown in Fig. 3, the two models produce very similar scores. In contrast, in Guangzhou, the LST shows minimal variation, while the DI exhibits noticeable fluctuations. The higher evaluation accuracy of DRSEI in Guangzhou (Fig. 3) suggests that DI differences may critically influence model performance. Notably, LST and DI spatial distributions correlate with model evaluation results: DRSEI and RSEI scores converge when their distributions align, whereas distribution divergences amplify score discrepancies between the models.

From Table 4, it is evident that the NDVI and the DI have relatively high loadings in both regions, corresponding to their significant spatial variations observed in Fig. 5. In contrast, the WET and the NDSI have lower loadings, consistent with their relatively minor spatial variations. Additionally, the LST exhibits greater variation in Beijing, leading to higher loadings. In contrast, Guangzhou’s LST shows less variation, resulting in lower loadings. This supports the reasonable assumption that an indicator’s loading magnitude depends on its variation extent within a region.

Although DI and LST maintain strong correlations with EI, he accuracy of DRSEI and RSEI differs significantly in Guangzhou. Specifically, the RSEI tends to overestimate EQ while the DRSEI underestimates it. This occurs because LST’s low variability in Guangzhou reduces its RSEI loading. Consequently, LST’s negative impact on EQ diminishes, elevating RSEI scores. Conversely, DI’s high variability in Guangzhou increases its DRSEI loading. This amplification of negative influence on EQ leads to lower DRSEI scores. This further demonstrates that the spatial heterogeneity of indicators directly affects their weighting in the models, ultimately influencing EQ evaluation results.

Spatial and temporal changes in ecosystem quality of the city proper in Chongqing

Given the superior evaluation accuracy of the DRSEI across various districts in Chongqing metropolitan area, this index was selected to a analyze changes in ecosystem quality within the region. To visually represent the spatial distribution changes in EQ, the DRSEI values were classified into five ecological EQ levels (Poor, Fair, Moderate, Good, and Excellent) at intervals of 0.2. The classification results are shown in Fig. 6.

As shown in Fig. 6, the ecological areas classified as ‘Moderate’ in 2002 largely transitioned to ‘Good’ by 2006. Between 2006 and 2010, regions with DRSEI values classified as ‘Moderate’ were primarily located in the western part of the study area, while a large portion in the eastern region improved to ‘Excellent’, indicating significant ecological enhancement in the east. However, from 2010 to 2014, the extent of ‘Excellent’ regions in the eastern part of the study area decreased substantially, with most of them downgraded to ‘Good,’ while the EQ in the central urban area further deteriorated to ‘Fair.’ This suggests that urban expansion may have exerted significant pressure on the ecological environment. Between 2014 and 2018, a large portion of the western area classified as ‘Good’ experienced slight ecological degradation, with most areas shifting from ‘Good’ to ‘Fair.’ From 2018 to 2022, the degradation trend persisted. The ‘Fair’ areas expanded westward, while the eastern region declined from ‘Good’ to ‘Moderate.’ These changes reflect an overall continuous decline in EQ.

To further quantify specific EQ changes in Chongqing’s city proper, the magnitude of changes was categorized into three levels (improved, basically stable, deteriorated) which were then expanded into five categories. The ‘improved’ category was further divided into ‘slightly improved’ and ‘markedly improved,’ while the deteriorated category was further divided into ‘mildly deteriorated’ and ‘markedly deteriorated.’ The specific classifications and their corresponding data were outlined in Table 5.

Table 5 shows that between 2002 and 2006, the ecological environment exhibited an overall improvement trend. The area classified as ‘improved’ significantly exceeded the ‘deteriorated’ area, with ‘slightly improved’ being the dominant category. This indicates an enhancement in EQ during this period. However, between 2006 and 2014, the trend reversed. The ‘deteriorated’ area markedly outpaced the ‘improved’ area, with ‘markedly deteriorated’ becoming the predominant classification. This reflects a sharp ecological decline during this period. From 2014 to 2022, although the ‘deteriorated’ remained larger than the ‘improved’ area, the gap between the two narrowed significantly compared to the previous period. Ecological degradation continued to be driven by ‘markedly deteriorated,’ while recovery primarily consisted of ‘slightly improved.’ Overall, 2002–2006 was the only period of overall EQ improvement in the study area, while 2006–2022 exhibited a continuous degradation trend. Notably, the extent of ecological deterioration was generally significant, whereas ecological recovery was relatively slow and limited. This suggests that over the past two decades, the recovery capacity of the ecosystem has struggled to offset the negative impacts of human activities and environmental changes.

To further reveal the spatial migration patterns of EQ, changes were visualized through ArcGIS software. Specific pathways of ecological pattern adjustments were explored, as shown in Fig. 7.

Overall, except for the ‘Excellent’ EQ level, which remained stably distributed in the southeastern Nanchuan District, the mean centers of all other levels were concentrated near the central area of the study region. Notably, the ‘Poor’ EQ level showed the strongest spatial volatility. It fluctuated near the study area’s center, with transient westward expansion to Rongchang District and eastward shift to Banan District. In contrast, other EQ levels exhibited constrained spatial displacements, mainly oscillating around the central area.

Although the spatial distribution of EQ levels underwent varying degrees of change between 2002 and 2022, by 2022 the mean centers of most levels had reverted to positions proximate to their 2002 baselines. This shift suggested two concurrent trends: the ecological pattern fluctuated over two decades while also exhibiting signs of stabilization. This may indicate that urban development and ecological management have prevented sustained expansion/contraction trends. Instead, spatial evolution adapted to existing ecological patterns.