Introduction

As the Industrial Revolution progresses, rapid economic expansion has increasingly been coupled with the degradation of ecological environments, posing a critical challenge for sustainable development. The question of how to balance economic growth with ecological preservation has become a pressing global concern1,2. By the end of the twenty-first century, the global urban area is expected to expand nearly sixfold, with Africa and Asia experiencing the most rapid urban growth3. Urban eco-environmental quality has gained significant attention in recent years due to the rapid urbanization processes occurring in both developed and developing countries4. The rapid pace of urbanization has led to various environmental challenges, such as unregulated land use change, agricultural land loss, air pollution, water pollution, and soil degradation5. Numerous studies have examined the spatiotemporal variations in urban eco-environmental quality, emphasizing the impact of urban expansion on local ecosystems6,7. Over the past four decades, China has emerged as the world’s fastest urbanizing nation, with its urban population steadily increasing8. In response to mounting ecological degradation and the need for sustainable economic development, the 19th National Congress of China has prioritized the construction of a “Beautiful China”9. The accelerated pace of urbanization has fostered a closer interaction between human activities and the environment, amplifying the negative impact of urbanization on ecological systems.

The development of remote sensing-based ecological evaluation models has played a pivotal role in monitoring and assessing the quality of the environment10. Initially, the Ministry of Environmental Protection of China introduced the Ecological Environment Index (EI) in 2006, which provided a numerical assessment of environmental quality11. However, EI was limited in its capacity to analyze the spatiotemporal variability of eco-environmental conditions. To address this gap, Xu Hanqiu introduced the RSEI in 2013, incorporating key factors such as soil, vegetation, surface temperature, and humidity12. The RSEI model, based solely on remote sensing data, offers an objective, robust, and comprehensive tool for environmental assessment13. Although the RSEI has proven to be an effective indicator for monitoring ecological quality at large scales, it has certain limitations, such as its sensitivity to land cover types and potential inaccuracies in areas with sparse vegetation. Additionally, RSEI may not fully capture the complexity of ecological processes in areas with heterogeneous environments, which could impact the precision of the results in certain regions. The limitation of RSEI should be emphasized in the process of ecological environment assessment.

Over the years, RSEI offers a valuable framework for assessing ecological quality across various regions, including urban environments14,15, mining areas16,17, and agricultural landscapes18. However, the application of RSEI in urban areas has gained particular attention due to the increasing urbanization worldwide. Recent advancements in RSEI methodology have incorporated additional factors such as aerosol optical depth (AOD) and socio-economic indicators, including population density, GDP, and land-use changes, to improve its predictive capabilities and better capture the complex interactions between urbanization and environmental quality19,20. These developments reflect the growing recognition that urban ecosystems are not solely shaped by natural processes, but are deeply influenced by socio-economic dynamics, policy decisions, and human activities21,22. Despite these improvements, few studies have investigated the spatiotemporal heterogeneity of eco-environmental quality across different urbanization stages, particularly during periods of rapid urban expansion. Understanding how urbanization stages impact ecological quality is critical for informed urban planning and sustainable development, especially in the context of emerging challenges such as urban heat islands, air pollution, and biodiversity loss.

The cities within the Beijing–Tianjin–Hebei–Shandong–Henan (BTHSH) region, which are at varying levels of urbanization, provide an ideal case study for this analysis. Understanding the spatiotemporal dynamics and driving forces of eco-environmental quality in these cities is crucial for addressing the structural contradictions between urban growth and environmental sustainability. Previous research has shown that both natural and anthropogenic factors, such as economic development23, urbanization24,25, and land use change26,27,28,29, significantly influence eco-environmental quality. Furthermore, the role of socio-economic factors in accelerating environmental degradation has become increasingly prominent30,31,32. The geographical detector method has emerged as a powerful tool for analyzing the interactions between multiple driving factors on eco-environmental quality33, revealing that the combined effects of these factors often have a more substantial impact than individual influences34,35. Geographical Detector to analyze the spatial heterogeneity of eco-environmental change, such as land use change36, ecosystem degradation37, and the effects of human activities on environmental health38.

This study aims to investigate the spatiotemporal heterogeneity and driving forces affecting eco-environmental quality across five cities in the BTHSH region using Landsat remote sensing data from 1985 to 2020, in conjunction with natural and anthropogenic datasets. This research explores the comparative spatiotemporal heterogeneity of ecological quality at different stages of urbanization in the Beijing–Tianjin–Hebei–Shandong–Henan region. This analysis provides novel insights into the impact of rapid urban expansion on ecological quality over time, a research area that has been largely overlooked in previous studies. By integrating long-term time series remote sensing data, socio-economic factors, spatiotemporal analysis, and predictions, our research offers a more comprehensive approach to assessing urban ecological quality, thereby laying the foundation for more informed urban planning and sustainable development strategies. Moreover, this study enriches the theoretical framework regarding the relationship between urban expansion and eco-environmental quality in rapidly urbanizing regions, a topic that has been relatively underexplored in the existing literature.

Materials and methods

Study area

Beijing, Tianjin, Hebei, Shandong and Henan (BTHSH) constitute the hub regions that connect northeast China, northwest China, southeast China, central China and overseas, between 31º23′ N ~ 42°37′ N, 110°21′ E ~ 122°43′ E (Fig. 1). It accommodates 21.6% of the country’s cities within 5.59% of the country’s land area. This makes it one of the largest urban-dense areas and the most densely populated area in China. Since the reform and opening up, according to the China Statistical Yearbook, China’s urbanization rate rose from 19.39% in 1980 to 53.73% in 2013 and is projected to reach 63.89% in 2020. By 2030, the urbanization rate is expected to range from 70 to 80%.The urbanization rate of the five cities in descending order is Beijing, Tianjin, Shandong, Hebei and Henan. Among them, the Beijing–Tianjin–Hebei urban agglomeration is China’s “capital economic circle”. Among them, Shandong Province has the largest coastal area followed by Hebei Province. Additionally, Since the establishment of the People’s Republic of China, urbanization has undergone three stages: the exploration and development stage from 1949 to 1979, the rapid development stage from 1979 to 2011, and the quality development stage from 2011 to the present (data collected through December 2024). To further analyze these stages, we identified key break points in the RSEI time series that correspond to significant shifts in urbanization. These break points are associated with major policy changes and urban development initiatives that coincide with the transitions between the phases of urban development. This analysis strengthens the relationship between urbanization stages and the evolution of the eco-environment. The selection of this study area is convenient to study the synchronization and correlation between urbanization rate and eco-environment evolution.

Fig. 1
figure 1

Geographic ___location and topographic features of the BTHSH study area.The map is generated using ArcGIS Pro 2.5 (ArcGIS Pro), URL: https://www.esriuk.com.

Full size image

Data source and processing methodology

At present, Landsat images are extensively utilized in eco-environment quality evolution and prediction research. This study utilized RSEI data (100-m resolution) provided by the National Earth System Science Data Center of China (https://gre.geodata.cn/). The dataset was generated based on the Google Earth Engine (GEE) cloud platform by collecting and processing Landsat TM/OLI/TIRS remote sensing imagery (30-m resolution) time series data from 1985 to 2020, with scenes selected to have less than 5% cloud cover in order to minimize atmospheric interference. The use of higher-resolution remote sensing data allows for a more accurate capture of subtle variations in land cover types, reducing uncertainties caused by land cover mixing effects. This is particularly beneficial in urbanized and complex geomorphological areas, where it enables a more refined identification of spatial variations in ecological quality. Remote sensing data of vegetation growth period were selected. In this study, RSEI data of the BTHSH regions in 1985, 1990, 1995, 2000, 2005, 2010, 2015 and 2020 were selected to analyze the evolution and driving mechanism of eco-environmental quality at different urbanization levels. The natural and urbanization-related data sources, which served as driving factors to analyze the underlying influences on RSEI changes, are presented in Table 1. These include DEM, NDVI, average annual rainfall, GDP, average annual temperature, population density (pop), and night light remote sensing datasets. The slope direction and slope data were extracted from DEM datasets. In the data preprocessing stage, the selected datasets projection was set to UTM-WGS84 coordinate system. In the process of geodetector data processing, all driver datasets were resampled to a spatial resolution of 1000 m using nearest neighbor resampling.

Table 1 Description of datasets used in this study, including spatial resolution and data sources.
Full size table

Methodology

RSEI eco-environment index

RSEI is constituted by four indicators, namely greenness, heat,wetness, and dryness, which represent ecological quality conditions, and all of them can be computed using remote sensing data. Greenness, heat, wetness, and dryness were represented by the vegetation index (NDVI), land surface temperature (LST), tasseled cap transformation moisture component (Wet), and the normalized difference building soil index (NDBSI), respectively, and their potential ecological associations were utilized to construct the RSEI (Fig. 2)35. The calculation of LST is based on Landsat thermal infrared band data, which is used to retrieve land surface temperature through a radiative transfer model. The calculation of NDBSI utilizes the reflectance difference between the Landsat red band and shortwave infrared band to assess the characteristics of buildings and bare soil.The principal component analysis (PCA) is applied in RSEI construction to synthesize NDVI, WET, NDBSI, and LST into a single index, ensuring objectivity and consistency without subjective weighting. To ensure comparability and eliminate unit inconsistencies, all component indicators (NDVI, LST, WET, and NDBSI) were rescaled to a uniform [0,1] range before being integrated via principal component analysis (PCA).The first principal component (PC1) is utilized to construct RSEI, and the contribution of each indicator to PC1 is regarded as the weight. Therefore, the RSEI index has the advantages of strong objectivity, robustness, accessibility, visualization, etc., and can be applied to various regions with different geographical conditions25,39. RSEI’s model can be briefly expressed as: RSEI = PC1[f(wetness, greenness, dryness, heat)]. The PC1 typically captures the majority of the variance and is used as the integrated ecological quality score.

Fig. 2
figure 2

RSEI construction.

Full size image

Spatial auto-correlation analysis

Spatial autocorrelation analysis is employed to assess whether the distribution of eco-environmental quality was clustered, including both global autocorrelation (global Moran’s I) and local autocorrelation(local Moran’s I). Global Moran’s I was employed to verify spatial element interrelationships in BTHSH region, with values ranging from − 1 to 1. When the global Moran’s I approaches 1, it suggests that the RSEI exhibits a significant positive correlation in space, that is, the numerical similarity of adjacent regions is high. When the global Moran’s I approaches − 1, it suggests a significant negative correlation among RSEI, namely, the values of adjacent regions are dissimilar. When the global Moran’s I approaches 0, it implies that the RSEI is not spatially correlated, that is, the distribution of the RSEI value among the regions occurs randomly. Local Moran’I is integrated into individual spatial units and local indicators of spatial association (LISA) for ecosystem quality are computed using a fishnet with a specific edge length40. Spatial autocorrelation analysis was conducted using ArcGIS Pro 2.5 software to evaluate the spatial clustering characteristics of RSEI values across the study area. The analysis was performed based on RSEI raster datasets with a spatial resolution of 100 m × 100 m.

$$Global \;Moran^{\prime}s I = \frac{{n\mathop \sum \nolimits_{i = 1}^{n} \mathop \sum \nolimits_{j = 1}^{n} W_{ij} \left( {x_{i} - \overline{x}} \right)\left( {x_{j} - \overline{x}} \right)}}{{\left( {\mathop \sum \nolimits_{i = 1}^{n} \mathop \sum \nolimits_{j = 1}^{n} W_{ij} } \right)\mathop \sum \nolimits_{i = 1}^{n} \left( {x_{i} - \overline{x}} \right)^{2} }}$$
(1)
$$\begin{array}{*{20}c} {Local \;Moran^{\prime}s I = \frac{{\left( {x_{i} - \overline{x}} \right)\mathop \sum \nolimits_{j = 1}^{n} W_{ij} \left( {x_{j} - \overline{x}} \right)}}{{\frac{1}{n} \times \mathop \sum \nolimits_{i = 1}^{n} \left( {x_{i} - \overline{x}} \right)^{2} }}} \\ \end{array}$$
(2)

In formula (1)–(2), W represents the space matrix; and \(\text{n denotes the observed value of the regional units}\), and \({\text{x}}_{\text{i}}\) is the observed value of the i-th unit. \(\overline{\text{x} }\) is the mean value of the observations. \({\text{W}}_{\text{ij}}\) is the spatial weight matrix for each cell.

The geographical detector

Geodetector is a methodology founded on mathematical statistics, which detects spatial heterogeneity in geographic phenomena and discloses the driving factors41. Due to the absence of homoscedasticity, normality and other assumptions, as well as the limitation of handling category variables, Geodetector has been extensively applied in meteorological, economic, ecological and other geography-related domains to identify the spatial heterogeneity of each influencing factor42,43,44. In this study, the factor and interaction detection modules of the Geodetector model were employed to identify the dominant factors influencing the spatiotemporal variation of RSEI and to assess how interactions between these factors affect changes in RSEI. The analysis was conducted using the Excel-based Geodetector tool. The factor detector quantifies the explanatory power of individual variables on the spatial distribution of RSEI, expressed as a q-statistic (0–1), with higher values indicating stronger influence. The interaction detector evaluates the combined effects of two variables, revealing whether their interaction enhances, diminishes, or independently affects spatial heterogeneity. This enables the identification of nonlinear synergistic effects between urban and environmental factors. In the interaction detector module of the Geographical Detector, interactions between pairs of factors are assessed based on their combined explanatory power regarding the dependent variable. If the joint explanatory power (q-value) of two factors is greater than that of each factor individually but less than the sum of their individual q-values, the interaction is classified as double-factor enhancement. If the joint q-value exceeds the sum of the individual q-values, the interaction is categorized as nonlinear enhancement, indicating a synergistic effect with more-than-additive influence on the spatial heterogeneity of eco-environmental quality.

The formula is:

$$\begin{array}{*{20}c} {q = 1 - \frac{{\mathop \sum \nolimits_{k = 1}^{L} N_{k} \sigma_{k}^{2} }}{{N\sigma^{2} }} = 1 - \frac{SSW}{{SST}}} \\ \end{array}$$
(3)
$$\begin{array}{*{20}c} {SSW = \mathop \sum \limits_{h = 1}^{L} N_{k} \sigma_{k}^{2} } \\ \end{array}$$
(4)
$$\begin{array}{*{20}c} {SST = N\sigma^{2} } \\ \end{array}$$
(5)

where h = 1, 2,…, n; L represents the partition of the independent variable and the dependent variable; N represents the total quantity of units; \({\text{N}}_{\text{k}}\) is the total quantity of units in layer h; \({\upsigma }_{\text{k}}^{2}\) is the variance of RSEI value corresponding to h layer.\({\upsigma }^{2}\) is the variance of RSEI value corresponding to the whole region. The SSW represents the aggregation of regional variances. The higher the q value is, the stronger the explanatory power of the corresponding factors becomes.

Methods of sustainability and anti-sustainability analysis

The Hurst index serves as an efficacious instrument for analyzing the long-term dependence of time series data. The long-term dependence refers to the continuous correlation of observed values in time series, which can help us understand the long-term trend of RSEI and predict the stability of the future trend of RESI. The Hurst index was initially proposed by the British hydrologist Hurst and has been extensively utilized in hydrology, geology, climatology, geology, economics, and other domains45. Based on the fundamental principle of the time series {RESI (t)}, For any positive integer τ ≥ 1, the mean sequence {\(\overline{{\text{RSEI}}_{(\uptau )}}\)} is defined as follows:

$$\overline{{RSEI_{\left( \tau \right)} }} = \frac{1}{\tau }\mathop \sum \limits_{t = 1}^{\tau } RSEI_{\tau } \quad \tau = 1,2 \ldots ,n$$
(6)

where t = 1, 2,… n. The deviation \({\text{X}}_{(\text{t},\uptau )}\), range \({\text{R}}_{(\uptau)}\) and standard deviation \({\text{S}}_{(\uptau )}\) of n observation points are respectively defined as:

$$X_{{\left( {t,\tau } \right)}} = \mathop \sum \limits_{t = 1}^{t} \left( {RSEI_{\left( t \right)} - \overline{{RSEI_{\left( \tau \right)} }} } \right) \quad 1 \le {\text{t}} \le \tau$$
(7)
$$R_{\tau = } max_{1 \le t \le \tau } X\left( {t,\tau } \right) - min_{1 \le t \le \tau } X\left( {{\text{t}},{\uptau }} \right) \quad \tau = 1,2, \ldots ,n$$
(8)
$$S_{\left( \tau \right)} = \left[ {\frac{1}{\tau }\mathop \sum \limits_{t = 1}^{\tau } (RSEI_{\left( t \right)} - RSEI_{\left( \tau \right)} )^{2} } \right]^{\frac{1}{2}} \quad \tau = 1,2, \ldots ,n$$
(9)
$$\begin{array}{*{20}c} {E\left[ {\frac{{R_{\left( \tau \right)} }}{s\left( \tau \right)}} \right] = A_{\tau }^{H} as \quad \tau \to \infty } \\ \end{array}$$
(10)

Using \(\text{s}\left(\uptau \right)\) to standardize \({\text{R}}_{\left(\uptau \right)}\), we get \({\text{R}}_{\left(\uptau \right)}\)/ \(\text{s}\left(\uptau \right)\); A is a constant; The Hurst index is denoted by H (H ∈ (0, 1)). When H is equal to 0.5, the RSEI sequence shows no autocorrelation. When 0.5 < H < 1, the RSEI sequence indicates a positive correlation, that is, sustainability. When 0 < H < 0.5, the RSEI sequence has a negative correlation, that is, anti-sustainability.

The Sen slope estimation method

The Sen slope estimation represents a robust non-parametric statistical trend calculation approach. This approach features high computational efficiency and is not sensitive to measurement errors as well as outlier data, and is appropriate for the trend analysis of long time series data of RSEI46. The Sen slope estimation regards the median slope as the overall trend of a time series. The formula is expressed as follows47,48:

$$\begin{array}{*{20}c} { \theta = Median\left( {\frac{{RSEI_{j} - RSEI_{i} }}{j - i}} \right) \quad \forall j > i} \\ \end{array}$$
(11)

where \({\text{RSEI}}_{\text{j}}\) and \({\text{RSEI}}_{\text{i}}\) are the values for year j and year i of the RSEI respectively. Where \(\uptheta\) is the slope which represents the variation trend of RSEI; Median() represents the median; if \(\uptheta\)> 0, it suggests an upward trend over time series; If \(\uptheta\) < 0, it implies a decreasing trend with time series, and the spatial change trend of RSEI in the BTHSH region can be analyzed accordingly.

Standard deviation ellipse method

The Standard Deviation Ellipse (SDE) method, originally introduced by D. Welty Lefever in 1926, is a classical spatial statistical technique used to analyze the directional characteristics and dispersion patterns of geographic phenomena49,50. SDE has become a routine statistical tool in ArcGIS Pro spatial statistics module, which can show the central trend, dispersion degree and direction trend of RSEI changes51. In this paper, the SDE tool of the ArcGIS Pro 2.5 software is applied to calculate two-dimensional RSEI data. Firstly, a new element class is created, and the ellipse is centered at the average center of all elements. The spatial distribution characteristics of RESI data center trend, dispersion and direction trend are represented by the long semi-axis, center point, short semi-axis, flatness and azimuth angle of the ellipse. The calculation formula for the aforementioned parameters is presented as follows:

$$M_{i} \left( {X_{i} ,Y_{i} } \right) = \left[ {\left( {{\raise0.7ex\hbox{${\mathop \sum \nolimits_{i = 1}^{n} W_{i} X_{i} }$} \!\mathord{\left/ {\vphantom {{\mathop \sum \nolimits_{i = 1}^{n} W_{i} X_{i} } {\mathop \sum \nolimits_{i = 1}^{n} W_{i} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\mathop \sum \nolimits_{i = 1}^{n} W_{i} }$}}} \right),\left( {{\raise0.7ex\hbox{${\mathop \sum \nolimits_{i = 1}^{n} W_{i} Y_{i} }$} \!\mathord{\left/ {\vphantom {{\mathop \sum \nolimits_{i = 1}^{n} W_{i} Y_{i} } {\mathop \sum \nolimits_{i = 1}^{n} W_{i} }}}\right.\kern-0pt} \!\lower0.7ex\hbox{${\mathop \sum \nolimits_{i = 1}^{n} W_{i} }$}}} \right)} \right]$$
(12)
$$tan\theta =\frac{(\sum_{i=1}^{n}{W}_{i}^{2}{\overline{X}}_{i}^{2}-\sum_{i=1}^{n}{W}_{i}^{2}{\overline{Y}}_{i}^{2})+\sqrt{{(\sum_{i=1}^{n}{W}_{i}^{2}{\overline{X}}_{i}^{2}-\sum_{i=1}^{n}{W}_{i}^{2}{\overline{Y}}_{i}^{2})}^{2}+4\sum_{i=1}^{n}{W}_{i}^{2}{\overline{X}}_{i}^{2}{\overline{Y}}_{i}^{2} }}{2\sum_{i=1}^{n}{W}_{i}^{2}{\overline{X}}_{i}^{2}{\overline{Y}}_{i}^{2}}$$
(13)
$${\sigma }_{x}=\sqrt{\frac{\sum_{i=1}^{n}{({W}_{i}{\overline{X}}_{i}\text{cos}\theta -{W}_{i}{\overline{Y}}_{i}\text{sin}\theta )}^{2}}{\sum_{i=1}^{n}{W}_{i}^{2}}}$$
(14)
$${\sigma }_{y}=\sqrt{\frac{\sum_{i=1}^{n}{({W}_{i}{\overline{X}}_{i}\text{sin}\theta -{W}_{i}{\overline{Y}}_{i}\text{cos}\theta )}^{2}}{\sum_{i=1}^{n}{W}_{i}^{2}}}$$
(15)

where \({\text{X}}_{\text{i}}\) and \({\text{Y}}_{\text{i}}\) are the central coordinates of each research unit; \({\text{W}}_{\text{i}}\) is the weight of the research unit. \({\text{M}}_{\text{i}}({\text{X}}_{\text{i}},{\text{Y}}_{\text{i}})\) represents the center of gravity coordinates of RSEI. θ is the elliptical orientation, that is, the angle formed by clockwise rotation from the true north direction to the major axis of the ellipse; while \({\upsigma }_{\text{x}}\) and \({\upsigma }_{\text{y}}\) are the standard deviations along the \(\text{X}\) and \(\text{Y}\) axes respectively.

Results and analysis

RSEI average trend analysis of time series

The trend analysis of the average RSEI from 1985 to 2020 (at 5-year intervals) as shown in Fig. 3 reveals distinct trends corresponding to key urban development phases. From 1985 to 2010, the RSEI values remained relatively stable, fluctuating between 0.517 and 0.522, suggesting a period of ecological stability. However, from 2010 onward, a noticeable decline in RSEI values occurred, dropping from 0.521 in 2010 to 0.496 in 2020, marking a significant inflection point.This inflection aligns with the transition in China’s urbanization process from rapid expansion (1979–2011) to a quality-oriented development phase (post-2011), as detailed in section “Study area”. Between 1985 and 2010, during the phase of rapid urbanization, eco-environmental quality in the BTHSH region demonstrated a mild upward trend, rising from 0.517 in 1985 to 0.521 in 2010. In contrast, during the period of urbanization improvement from 2010 to 2020, the RSEI declined, indicating the complexity of improving ecological quality after rapid urbanization.

Fig. 3
figure 3

RSEI average trend analysis (1985–2020).

Full size image

The results highlight a significant turning point around 2010, correlating with policy shifts and socio-economic changes in China’s urban development trajectory. During the phase of rapid urbanization, the urbanization level advanced rapidly, resulting in RSEI values entering a lower range. However, after transitioning to the stage of urbanization improvement, the RSEI values did not exhibit immediate significant improvement, indicating that the recovery and enhancement of ecological environmental quality require a certain transitional period. The changes in ecological environmental quality exhibit a certain lag effect, meaning that environmental improvement measures or policy adjustments in the urbanization process typically take time to reflect in indicators such as RSEI. This lag effect indicates that, despite the implementation of improvement measures in certain stages, the recovery of ecological environmental quality does not occur immediately but requires a period of time to accumulate and manifest its effects.

Spatiotemporal variation characteristics of ecosystem evaluation based on RSEI

Based on eight RSEI images from 1985, 1990, 1995, 2000, 2005, 2010, 2015, and 2020, the RSEI of Beijing, Tianjin, Hebei, Shandong, and Henan province (municipality) was classified at an interval of 0.2. The reclassification of RSEI was divided into five indicators: excellent (0.8–1], good (0.6–0.8], moderate (0.4–0.6], fair (0.2–0.4], and poor (0, 0.2]. The spatiotemporal distribution of the RSEI in BTHSH regions from 1985 to 2020 reveals significant regional and temporal variations in ecosystem quality. Since the 1980s, China’s eco-environment has also undergone remarkable changes under the impact of comprehensive factors such as the rapid urbanization. It can be observed from Fig. 4 that in the BTSHS regions during the period from 1985 to 2020, the quality of the eco-environment is deteriorating on the whole as the urbanization level increases. In the BTSHS regions from 1985 to 2020, the proportions of the areas with a moderate level of eco-environment assessment and those with good and excellent levels demonstrated a decreasing tendency from 90.44% to 83.06%, while the area with poor and fair levels demonstrated an increasing trend from 9.56% to 16.94%. However, the percentage of the category area for the moderate type ecological environment index exhibited negligible variation, with the maximum value reaching 78.38%, the minimum value being 68.65%, and the average value standing at 73.42%. It can also be observed from Fig. 5 that there is a trend of deterioration in the eco-environment quality. According to Fig. 4, the RSEI of the BTHSH regions during 1985–2020 is mainly in the moderate classes, and the good and excellent ranges are mainly distributed in coastal cities. However, after 2010, the area of coastal cities with good eco-environment quality witnessed a significant decline and turned into the moderate classes. This is consistent with the research conclusion that from 1985 to 2020, the economic development of coastal cities in China is more rapid than that of inland regions, and the habitat quality of coastal cities is significantly degraded under the influence of rapid urbanization and industrialization52,53.

Fig. 4
figure 4

the percentage of areas classified by different eco-environment quality in the BTHSH region.

Full size image
Fig. 5
figure 5

Spatial distribution of eco-environmental quality classification in BTHSH region from 1985 to 2020.The map is generated using ArcGIS Pro 2.5 (ArcGIS Pro), URL: https://www.esriuk.com.

Full size image

As can be observed from Figs. 6, 7, 8, 9 and 10, the maximum reduction in the percentage of the eco-environment area at the moderate level in Beijing, Tianjin, Shandong, Hebei and Henan is respectively 29.3%, 20.91%, 8.72%, 16.79% and 14.34%, which is basically similar to the ranking of urbanization levels in each region. Although the area percentage of moderate level eco-environment quality in Shandong Province especially the eastern coastal cities (Fig. 5) decreased slightly from 1985 to 2020, the area percentage of the good classes dropped from 31.93% to 15.99%, and the area percentage of the poor classes rose from 4.41% to 17.37% (Fig. 9). With the advancement of urbanization, the land cover of coastal cities changes dramatically, and it has been demonstrated that urbanization has a more serious impact on coastal cities54,55. The maximum increases of the eco-environment assessment of the fair and poor classes in Beijing, Tianjin, Shandong, Hebei and Henan were 16.42%, 21.55%, 13.4%,8.54% and 8.09%, respectively, which is basically similar to the ranking of urbanization levels in each region. Beijing, as the political and economic center of the country, has a lower degree of eco-environment deterioration compared with Tianjin, followed by Shandong, Hebei and Henan. As the capital, Beijing is characterized as the country’s political center, cultural center, international communication center, and scientific and technological innovation center. According to the “Beijing Urban Master Plan (2016–2035)”, the objective is to initially establish a world-class harmonious and livable city by 203556. Tianjin serves as the national advanced manufacturing research and development base, the core area of northern international shipping, the financial innovation operation demonstration zone, and the reform pilot demonstration zone. Its role in financial innovation and advanced manufacturing has been prominently highlighted in the Outline of the Plan for the Coordinated Development of the Beijing–Tianjin–Hebei Region57. Therefore, although the urbanization level of Beijing is superior to that of Tianjin, the proportion of the area where the eco-environment quality has deteriorated in Beijing is marginally lower than that of Tianjin.

Fig. 6
figure 6

The percentage of areas classified by different eco-environment quality in the Beijing region.

Full size image
Fig. 7
figure 7

The percentage of areas classified by different eco-environment quality in the Tianjin region.

Full size image
Fig. 8
figure 8

The percentage of areas classified by different eco-environment quality in the Hebei region.

Full size image
Fig. 9
figure 9

The percentage of areas classified by different eco-environment quality in the Shandong region.

Full size image
Fig. 10
figure 10

The percentage of areas classified by different eco-environment quality in the Henan region.

Full size image

The changing trend of the RSEI standard ellipse

To investigate the spatial distribution trend of the eco-environment in the BTHSH regions, the standard deviation ellipse of RSEI for eight years (1985, 1990, 1995, 2000, 2005, 2010, 2015, and 2020) was constructed. The standard elliptic analysis approach was employed to explore the spatial distribution trend of eco-environmental quality in-depth. As can be observed from Fig. 11, the dominant direction of spatial distribution of eco-environment quality has little change, and the long axis of the ellipse shows a northeast-southwest direction, indicating that the spatial distribution pattern of ecological environment is northeast-southwest. From 1985 to 2020, the dominant direction of the eco-environment distribution trend in the BTHSH study area did not altered significantly, but its eco-environment center of gravity, agglomeration and diffusion characteristics gradually changed with time (Fig. 11). Firstly, the eco-environment elliptic center of the study area is situated in Shandong Province, adjacent to Hebei and Henan provinces, and the elliptic center demonstrates a trend of shifting from southwest to northeast towards the Beijing–Tianjin–Hebei capital economic circle. It suggests that the center of gravity of the eco-environment is gradually moving towards the area with a higher urbanization rate. The elliptic oblate demonstrated a tendency of increasing fluctuation, and the elliptic oblate increased significantly from 0.41 in 1985 to 0.58 in 2020 (Table 2). The increase in oblateness also suggests that the spatial distribution of eco-environment quality becomes increasingly directional with the escalation of the urbanization rate from 1985 to 2020. This law can also be found from the change of ellipse area. The area of the ellipse remained relatively unchanged, and the distribution of eco-environment quality presented relatively diffuse spatial characteristics (Fig. 11). Nevertheless, from 2015 to 2020, the area of the ellipse decreased significantly, and the distribution of eco-environment quality exhibited aggregation and enhanced directivity.

Fig. 11
figure 11

Standard ellipse variation trend of eco-environmental quality in the BTHSH regions.The map is generated using ArcGIS Pro 2.5 (ArcGIS Pro), URL: https://www.esriuk.com.

Full size image
Table 2 Average value and standard elliptical parameters of eco-environmental quality in the BTHSH regions.
Full size table

By integrating the average value of RSEI and the standard ellipse, the influences of different urbanization stages on the quality of the eco-environment were analyzed. Based on the average RSEI value of the BTHSH regions from 1985 to 2020, the overall eco-environment quality of the study area demonstrated a trend of slight fluctuation and decline over the past 35 years. Among them, the maximum value was 0.522 in 2000, and the minimum value was 0.496 in 2020 (Table 2). However, in the improvement stage of urbanization, the RSEI center of gravity of eco-environment quality increased from 0.408 in 2010 to 0.661 in 2020, and the ellipticity ascended from 0.38 to 0.58 (Table 2). During the second stage of urbanization development, the average value of eco-environment quality slightly declined, yet the spatial distribution trend became more pronounced, and the area with a superior ecological environment shifted towards the Beijing–Tianjin–Hebei direction.

Spatial autocorrelation analysis on the variations of eco-environmental quality

To further explore the spatial correlation and differentiation features of eco-environmental quality, both the global Moran’I and local Moran’I were employed to analyze the eco-environmental quality in the BTHSH regions. From 1985 to 2020, the global Moran’I was consistently greater than 0 (P ≤ 0.01 in the significance test), indicating that RSEI in the BTHSH regions exhibited a spatial positive correlation and presented a significant clustering feature in space (Table 3). During the period of rapid urbanization from 1985 to 2010, the global Moran’I of the BTSHS regions exhibited an ascending trend, ranging from 0.729 in 1985 to 0.765 in 2010. The improvement stage of urbanization from 2010 to 2020 witnessed a decreasing trend of the global Moran’I index, dropping from 0.729 in 2010 to 0.318 in 2020.In the second stage of urbanization quality improvement, the urbanization rate increases slowly, and the positive correlation of eco-environment quality is weakened by the interference of ecological protection policies and human activities, mitigating the adverse impacts of urbanization on the quality of the eco-environment.

Table 3 The global Moran’I of eco-environment quality in the BTHSH region.
Full size table

As depicted in Fig. 12, the local spatial autocorrelation of RSEI in BTHSH regions from 1985 to 2020 can be employed to analyze the specific locations of aggregation or dispersion features. The rapid urbanization also exerted an impact on the local spatial autocorrelation of the eco-environment. From 1985 to 2020, the High-High cluster types are predominantly distributed in the eastern coastal cities of Shandong Province, while the Low-Low cluster areas are mostly located in Henan Province and the southwestern part of Shandong Province. However, Beijing and Tianjin, which have a relatively high level of urbanization, did not have a local autocorrelation aggregation area of eco-environment quality. After 2000, the local spatial cluster characteristics of the Beijing–Tianjin–Hebei capital economic circle showed no significant characteristics, which was consistent with the conclusion that the global spatial autocorrelation was weakened. Meanwhile, the Low-Low cluster area in Hebei Province disappeared after 2000 as well. Additionally, in line with the advancement of urbanization, the concentrated area of eco-environment quality demonstrates a tendency of fluctuating in a decreasing manner. During the period of urbanization improvement from 2010 to 2020, the area of High-High cluster areas of coastal cities in Shandong Province has decreased significantly, which is consistent with the conclusion that the area of coastal cities with more serious urbanization interference and good classes of eco-environment quality has decreased significantly.

Fig. 12
figure 12

Local spatial autocorrelation distribution of BTSHS region from 1985 to 2020. The map is generated using ArcGIS Pro 2.5 (ArcGIS Pro), URL: https://www.esriuk.com.

Full size image

Forecast of the changing trend of eco-environment quality and its driving factors in the BTHSH region

Theil-sen median trend degree of RSEI from1985 to 2020

By analyzing the variation of the Theil-Sen Median trend degree of RSEI (Fig. 13), it can be observed that from 1985 to 2020, the area of RSEI with sen > 0 accounted for 46.51% in the BTHSH region, suggesting that the eco-environment in this region demonstrated a favorable trend over the time series, and the area with sen < 0 accounted for 53.88%. The results indicated that from 1985 to 2020, along with the rapid advancement of urbanization, the overall eco-environment in this area showed a downward trend, particularly in coastal cities. From the aspect of spatial distribution (Fig. 13), the eco-environment of the ecological conservation area in the northwest of Beijing–Tianjin–Hebei demonstrated a distinct tendency of improvement, while the eco-environment quality of the eastern coastal development area presented an obvious trend of deterioration. In the Beijing–Tianjin–Hebei capital economic circle, the proportion of the area with sen < 0 is 20.59%, and that of the area with sen > 0 is 19.65%. The overall RSEI of Henan and Shandong province also showed a small downward trend, with 33.3% of the area with sen < 0 and 26.48% of the area with sen > 0.The eco-environment in the eastern coastal area of Shandong province showed an obvious downward trend (Fig. 13). The overall eco-environmental quality of Henan Province demonstrates a favorable trend, with sen < 0 area accounting for 14.2% and sen > 0 area accounting for 16.91%, which is closely related to the fact that Henan’s cultivated land area accounts for 64.5%, which is the province with the largest proportion of cultivated.

Fig. 13
figure 13

Theil-Sen Median trend distribution of RSEI from1985 to 2020.The map is generated using ArcGIS Pro 2.5 (ArcGIS Pro), URL: https://www.esriuk.com.

Full size image

Analysis of the future changing trend of eco-environment quality

Hurst index is used to determine the persistence or anti-persistence of a trend component. When the index H equals 0.5, it implies that RSEI exhibits no temporal autocorrelation and possesses randomness. When 0 < H < 0.5, it suggests that there is a negative temporal correlation, showing anti-persistence. When 0.5 < H < 1, it suggests that there is a positive temporal correlation and a persistence. As is clearly demonstrated in Fig. 14a, most regions of the BTHSH have a persistence in time series. The area where the Hurst index is greater than 0.5 constitutes 84.1%, while the area where the Hurst index is less than 0.5 amounts to 15.9%.

Fig. 14
figure 14

Hurst index and the future trend variation of RSEI. (a) The RSEI continuity and anti-continuity distribution; (b) the future trend variation of RSEI.The map is generated using ArcGIS Pro 2.5 (ArcGIS Pro), URL: https://www.esriuk.com.

Full size image

Figure 14b presents the future trend variation of RSEI obtained through the superimposition analysis of the Theil-Sen Median trend degree change and the Hurst index. When sen > 0 and hurst > 0.5 or sen < 0 and hurst < 0.5, RSEI demonstrated a favorable trend. When sen < 0 and hurst > 0.5 or sen > 0 and hurst < 0.5, RSEI exhibits a variable trend.Therefore, the hurst index is reclassified into two categories with 0.5 as the middle value. Meanwhile, the value range of the sen trend degree [− 0.053, 0.095] is evaluated. Based on the natural breakpoint method and the extent of increase or decrease of the RSEI, it is reclassified into seven categories. The grid multiplication operation of these two reclassification data is employed to acquire seven distinct trends of RSEI in the future. Based on this, the future change trend of RSEI can be classified into seven types: significant variation, moderate variation, slight variation, basically unchanged, slight improvement, moderate improvement, and obvious improvement.

It can be observed from the Table 4 and Fig. 14b that regarding the future changing trend of the BTHSH regions, the area of basically unchanged types holds the highest proportion, suggesting that the overall trend of the RSEI index will not undergo significant alterations in the future. The proportion of the area of the basically unchanged type in the five provinces (municipalities) in the future is 69.82%, 63.70%, 66.29%, 63.83% and 62.66%, respectively. Among the slightly improved types, the area of Henan Province takes up the highest proportion of 23.58%, followed by Hebei province with 21.32%. The coastal city of Tianjin and Shandong province have relatively low proportions of the area, accounting for 12.41% and 12.81% respectively, and Beijing, as the capital city, accounts for 16.14%. As can be observed from the spatial distribution map (Fig. 14b), the slight improvement trend is mostly distributed in the central and southeast parts of Hebei, the north and southeast parts of Henan and the southwest and north parts of Beijing. The types of slight variation of the future trend were mostly distributed in the southeast of the coastal city of Tianjin and the eastern parts of Shandong Province, accounting for 22.87% and 22.16% respectively. The area proportion of the slight variation type in Beijing and Henan province is 12.17% and 12.81% respectively, while the area proportion of the slight variation type in Hebei Province is the lowest, being 11.15%. There are few types of moderate change and moderate improvement in the future, among which Beijing has the largest proportion of moderate variation (1.5%) and Hebei province has the largest proportion of moderate variation (0.77%). The types of obvious improvement and significant variation in the future are merely distributed in the capital city of Beijing, accounting for 0.01% and 0.03% respectively, which have a significant correlation with the interference intensity of human activities.

Table 4 RSEI future change trend area percentage in BTHSH regions.
Full size table

Driving mechanism of eco-environmental quality change in the BTHSH regions

This study comprehensively considers the availability and scientificity of factors and selects GDP, light night data, and population density data as the factors reflecting urbanization58, as well as temperature, rainfall, NDVI, elevation, slope, and slope direction as the natural factors59,60. Night light data, as an important means to characterize the urbanization process and human activities, has been proved to be closely related to the economic scale, urban population and expansion61.

According to the single factor analysis of geographical detector, nine independent variables of natural and human factors, such as NDVI, temperature and GDP, have substantial impacts on the spatial differentiation of the dependent variable RSEI. By analyzing the accumulation bar chart of q values (Fig. 15) for the single factor detection results from 1985 to 2010, it is discovered that the driving factors such as temperature and GDP rank first to second, followed by rainfall, DEM, NDVI, population density, night light index, slope and slope direction(Significance test P < 0.05)62. During The second stage of urbanization development from 2010 to 2020, NDVI, night light index, POP, and GDP were ranked among the top four driving forces successively, and the contribution of social and urbanization related data increased in comparison with that of natural driving factors.

Fig. 15
figure 15

The accumulation bar chart of q values for the single factor detection from1985 to 2020.

Full size image

As is shown in Fig. 15, during the first stage of urbanization development from 1985 to 2010, the q values of GDP, NDVI and population density (POP) all demonstrated a downward trend. However, in the second stage, they exhibit an increasing trend. The driving factors of light night index, which represent the level of urbanization, present a fluctuating upward trend throughout the entire stage of urbanization development, with the q value rising from 0.019 in 1985 to 0.049 in 2020. While natural factors driving force contribution temperature and rainfall all demonstrate a trend of decreasing fluctuation63. The q value of temperature dropped from 0.059 in 1985 to 0.005 in 2020, and the q value of rainfall declined from 0.066 in 1985 to 0.017 in 2020. The influence of DEM, slope and slope direction on the eco-environmental quality of the BTHSH regions has little change. The impact of human factors such as the light night index and population density demonstrated a fluctuating growth trend from 1985 to 2020, suggesting that human factors such as social economy and urbanization have increasingly significant impacts on eco-environment64. Similar conclusions were reached through the study of Taiyuan65, a coastal city in China66, the Qinghai-Tibet Plateau58, and other regions. Their research indicates that natural factors serve as the primary elements affecting the quality of the eco-environment, whereas the influence of human activities such as urbanization and economic activities is growing increasingly prominent.

To sum up, during the rapid urbanization process in China from 1985 to 2020, human factors such as GDP and the night light index have increasingly become significant drivers of the alterations in the spatial distribution pattern of RSEI particularly in the second stage of urbanization67. During the first phase of urbanization spanning from 1985 to 2010, natural factors like temperature and precipitation were the principal elements influencing RSEI. In the second stage of urbanization development from 2010 to 2020, driven by the concept of sustainable development in China, the influence of human disturbance factors, such as the NDVI, economic factors, population density, and the light night index on RSEI has gradually increases68.

Investigating the interaction factors driving the evolution of the eco-environment

The driving force of the temporal and spatial change of RSEI was analyzed by the geographical detector factor interaction (Fig. 16). The interaction among influencing factors manifested in two primary types: double-factor enhancement and nonlinear enhancement. This indicates that the combined effects of certain factor pairs were more pronounced than their individual contributions. Consequently, it can be inferred that changes in eco-environmental quality in the BTHSH region are not solely attributable to isolated factors but rather stem from intricate interactions, particularly synergistic relationships, among topography, climate, vegetation, and socio-economic elements. As depicted in Fig. 16, the factors corresponding to the independent variables are as follows: X1 is DEM, X2 is GDP, X3 is NDVI, X4 is slope, X5 is population density(POP), X6 is slope direction, X7 is the annual average temperature, X8 is the annual precipitation, and X9 is the night light index. From 1985 to 1990, the interaction among temperature, rainfall, and other factors had a significant impact on the eco-environment, especially the interaction between temperature and GDP, which ranged from 0.21 to 0.15. During the period from 1995 to 2010, the interaction between temperature and other factors exerted a considerable influence on the eco-environment, especially the interaction between temperature and light night index changing from 0.06 to 0.18. From 1985 to 2020, particularly after 2000, the interaction between the night light index and other data has been strengthened, suggesting the growing influence of urbanization on the eco-environmental quality. Since 2015, the influence of multi-factor interaction on the eco-environment has been growing, and the contribution of comprehensive factor interaction to the alteration of the eco-environment has increased. Furthermore, the connection between the eco-environment and the influencing factors has become more complex69. During the first phase of urbanization development from 1985 to 2010, the interaction q value of annual mean temperature and annual precipitation was relatively high, ranging from 0.19 to 0.26, signifying that the interaction between the two exerted a significant impact on the eco-environmental quality in the study area70. During the second phase of urbanization development from 2010 to 2020, the q value of interaction between night light index and other factors gradually increased from 0.03 to 0.10.The interaction q value of GDP, population density(POP) also demonstrated an upward trend in the second phase of urbanization development. The interaction effect of night light index with other natural and economic factors gradually strengthened, indicating that the implementation of ecological protection policies and urbanization development have an increasingly significant influence on the eco-environmental quality in the study area71,72.

Fig. 16
figure 16

The q value of RSEI evolution driving factors interaction detection in BTHSH regions from 1985 to 2020.

Full size image

Discussion

Comparison with previous studies

Currently, there is a lack of long-term quantitative comparative studies on the spatiotemporal evolution of eco-environmental quality across cities at various urbanization stages in the BTHSH region, particularly focusing on the capital economic circle and its surrounding cities. This study adopts a long-term temporal analysis (1985–2020) with five-year intervals to assess ecological quality dynamics in the BTHSH region73,74. Unlike previous studies relying on limited or coarse time points, this approach enhances temporal sensitivity, enabling detection of both immediate responses and cumulative effects of urbanization. The inclusion of the 1985–2000 period further strengthens the analysis by capturing early ecological trends, offering a more robust understanding of spatiotemporal evolution across different urbanization stages75,76,77. This extended and refined temporal framework provides novel insights into the long-term impacts of rapid urban growth on ecological quality, highlighting the importance of continuous monitoring and fine-scale temporal resolution in urban ecological research78. The observed variations in eco-environmental quality across the BTHSH region highlight that ecological change is rarely driven by a single factor in isolation. Instead, it reflects the result of complex, often synergistic, interactions among natural and socio-economic variables. These include the combined influences of topographic features, climatic conditions, vegetation dynamics, and human activities79,80. Moreover, the influence of urbanization and the implementation of environmental protection policies on RSEI exhibits a notable lag effect, suggesting that improvements or degradation in ecological quality often do not occur immediately, but instead emerge gradually over time following changes in urban development and governance78.

Limitations and future research direction

The RSEI offers an efficient and scalable methodology for ecological quality assessment, yet it is subject to certain limitations. Compared with conventional ecological indicators, such as the NDVI, vegetation coverage, or composite indices based on subjectively assigned weights, the RSEI provides a more integrated and objective method by synthesizing four critical ecological dimensions (greenness, wetness, dryness, and heat) using PCA81,82. These indicators are directly extracted from satellite imagery and can be reliably obtained across extensive spatial domains and extended temporal scales, rendering RSEI particularly appropriate for regional-scale dynamic assessments. Although the RSEI is a widely utilized composite index for ecological quality assessment, it reflects a relative ecological status and may not comprehensively capture all local environmental characteristics, especially in regions with complex ecosystems or substantial human disturbances83,84,85. Furthermore, the construction of RSEI depends on PCA, which, despite being an objective and data-driven approach, possesses inherent limitations. Future research could delve into the integration of auxiliary data, including land use patterns, socio-economic factors, and biodiversity indicators, while exploring advanced data fusion methodologies to enhance the comprehensiveness and adaptability of the ecological environment quality assessment framework.

This study identifies a complex and dynamic relationship between eco-environmental quality and its driving factors in the BTHSH region, where natural factors traditionally play a dominant role, while urbanization increasingly influences ecological changes. However, the explanatory power of individual factors fluctuates across different regions and time periods, indicating the inherent complexity of these interactions. For instance, in the first batch of civilized pilot areas in Jiangxi Province, factors such as slope, population density (POP), and elevation are critical to ecological quality14, while studies on mining cities in China highlight the importance of factors such as mining distance, digital elevation model (DEM), and GDP in shaping eco-environmental outcomes86. Such variability emphasizes that the drivers of eco-environmental change are highly context-specific, varying across different regions and time periods. Specifically, although the use of high-resolution (30-m) remote sensing data has improved the spatial accuracy of ecological quality assessments, certain limitations remain in capturing fine-scale ecological heterogeneity, particularly in highly urbanized and transitional landscapes. Future studies could further refine ecological assessments by integrating higher-resolution remote sensing data to better delineate localized environmental variations.

Additionally, over the past two decades, significant ecological initiatives, including the "Regional Major Strategic Ecological Environment Protection Strategy"87 and the "National Ecological Environment Protection Outline"88, have accelerated environmental protection efforts and mitigated some of the negative impacts of urbanization. These policy interventions have been effective in enhancing ecological quality, as evidenced by the increase in areas with “good” eco-environmental quality between 2015 and 2020. Nevertheless, the direct impact of these policies on the quantitative aspects of eco-environmental quality remains difficult to isolate within the driving force model, due to the complex interplay of multiple factors and the time lag between policy implementation and observable environmental changes. It is the challenge of fully capturing the influence of all potential drivers, particularly non-quantitative or less easily measurable factors such as policy changes, social behavior, and institutional factors82,89,90,91,92,93. Although this study incorporates a range of natural, economic, and social variables, future research should aim to integrate a more comprehensive set of socio-political factors and explore the causal pathways in greater depth. As research progresses, the inclusion of additional influential factors and advanced analytical techniques will enable a more nuanced understanding of the eco-environmental dynamics in the BTHSH region and beyond. This could offer valuable insights into optimizing environmental governance and formulating more targeted strategies for sustainable urbanization.

Conclusion

This study employed the RSEI, constructed using GEE and long-term remote sensing data, to assess the eco-environmental quality of the BTHSH (Beijing–Tianjin–Hebei–Shandong–Henan) region from 1985 to 2020 across different stages of urbanization. By utilizing tools such as standard ellipse, spatial autocorrelation analysis, the Hurst index, and the Theil-Sen Median trend, we analyzed both the historical and future spatiotemporal evolution of eco-environmental quality, as well as the driving mechanisms underlying these changes. The main findings of the study are as follows:

  1. (1)

    Spatiotemporal Heterogeneity of Eco-environmental Quality.

    The results indicate a general deterioration in eco-environmental quality as urbanization progresses. From 1985 to 2010, during the early phase of urbanization, the eco-environmental quality showed a fluctuating but overall upward trend, increasing from 0.517 to 0.521. However, in the second urbanization phase (2010–2020), eco-environmental quality slightly declined, with a more pronounced shift in spatial distribution, particularly towards the Beijing–Tianjin–Hebei region. The Global Moran’s I statistic decreased from 0.729 in 1985 to 0.318 in 2020, suggesting a weakening of spatial autocorrelation as urbanization intensified. Notably, coastal cities, especially in Shandong Province, experienced a significant deterioration, with areas of poor eco-environmental quality increasing from 8.15% to 17.94% during this period.

  2. (2)

    Future prediction on the quality of eco-environment.

    Forecasts indicate that Henan Province will see the most significant improvement in eco-environmental quality (23.58%), followed by Hebei (21.32%). In contrast, coastal cities like Tianjin and Shandong are projected to experience relatively smaller improvements (12.41% and 12.81%, respectively), while Beijing, as the capital, is expected to see a 16.14% improvement. This highlights a pronounced negative impact of urbanization on the eco-environmental quality of coastal cities, emphasizing the need for targeted environmental management strategies in these areas.

  3. (3)

    Driving forces of eco-environmental quality evolution.

    The analysis of driving factors reveals that, during the early urbanization phase (1985–2010), temperature, GDP, and rainfall were the primary drivers of eco-environmental quality, with temperature having the greatest explanatory power. In contrast, during the later urbanization phase (2010–2020), indicators such as NDVI, the night light index, population density (POP), and GDP emerged as the most significant factors. The increasing interaction between the night light index and other socio-economic and environmental factors underscores the growing influence of urbanization on eco-environmental quality in the study region.

This research offers critical insights into the spatiotemporal dynamics of eco-environmental quality in rapidly urbanizing regions, providing valuable information for policymakers and urban planners. The findings underscore the importance of considering both natural and socio-economic factors in sustainable urban development strategies. The study also highlights the need for region-specific approaches to mitigate the adverse effects of urbanization, particularly in coastal cities, to ensure the long-term sustainability of urban ecosystems.