Multiscale dynamics and nonlinear associations of ecological quality in the Hangzhou Bay urban agglomeration, China

Ecological quality change in rapidly urbanizing regions is spatially heterogeneous, but assessments based on a single analytical scale may obscure local changes and alter the apparent factor importance. This study constructed a 2000–2024 Remote Sensing Ecological Index (RSEI) series for the Hangzhou Bay Urban Agglomeration using Google Earth Engine and Landsat imagery. Ecological changes were compared across 1, 3, 5, and 10 km grids and county units with XGBoost–SHAP model fitted to pooled observations from six benchmark years. Mean RSEI showed an early decline followed by sustained recovery, yielding a net increase of 0.0423. High values were concentrated in southwestern forested areas, whereas low values occurred mainly in the northern plain, densely urbanized bay areas, and urban-expansion corridors. The proportion of improved area increased from 75.34% at 1 km to 84.25% at the county scale. During 2015–2020, county-level aggregation increased the stable proportion while reducing the improved proportion, indicating that coarse units can mask opposing local changes. Land-use intensity had the largest relative mean absolute SHAP contribution at all scales (44.59%–55.43%), while the contribution of nighttime light intensity rose with aggregation and reached 20.33% at the county scale. These results demonstrate that ecological-change detection and interpretation of modeled associations are scale dependent, supporting a nested multiscale assessment framework.

Rapid urbanization is profoundly reshaping regional land-use patterns, the surface energy balance, and ecological processes, contributing to habitat fragmentation, an intensified thermal environment, reduced hydrological regulation capacity, and mismatches between the supply of and demand for ecosystem services (Liang et al., 2025; Zhao et al., 2025; Zhu et al., 2025). These pressures are especially complex in urban agglomerations, where intensive urban construction, population and industrial concentration, transport-network expansion, and the reconfiguration of blue–green spaces interact to produce pronounced spatial heterogeneity, stage-specific fluctuations, and scale dependence in ecological quality (Liu et al., 2025; Lu et al., 2025). Ecological quality is not only an important indicator of regional sustainability but also a basis for identifying ecological degradation risks, evaluating restoration effectiveness, and improving territorial spatial governance. A continuous, objective, and comparable monitoring framework is therefore essential for characterizing ecosystem changes in rapidly urbanizing regions and examining their associated factors.

Remote sensing provides essential support for large-area, long-term monitoring of ecological quality (Lima et al., 2025; Wulder et al., 2022). Compared with individual vegetation indices or land-use indicators, the remote sensing ecological index (RSEI) integrates greenness, wetness, dryness, and heat and uses principal component analysis to reduce subjectivity in assigning weights, thereby providing a more comprehensive representation of surface ecological conditions (Xu et al., 2019). RSEI has been applied to ecological monitoring in cities, watersheds, mining areas, protected areas, and urban agglomerations (Sun et al., 2022; Wu et al., 2025; Yuan et al., 2021). With the development of cloud-computing platforms such as Google Earth Engine (GEE), constructing long-term RSEI series from multisource remote sensing data, including Landsat imagery, has become an important research direction (Gorelick et al., 2017). Previous studies have found that regional ecological improvement can coexist with local degradation in core built-up areas and along urban expansion fronts in rapidly urbanizing regions (Guo et al., 2024; Zhang et al., 2022). Nevertheless, most studies remain focused on describing spatiotemporal patterns, with limited attention to cross-year comparability in long-term series, the consistency of change detection across spatial scales, and the scale sensitivity of model-derived factor contributions and interactions (Guo et al., 2025; He et al., 2025a).

Spatial scale is a central consideration in both identifying ecological-quality patterns and interpreting associated factors (Comber and Harris, 2022; Liu et al., 2023a). Vegetation growth, topographic constraints, and hydrothermal conditions often operate over relatively broad areas, whereas road proximity, urban expansion, local thermal conditions, and land-use disturbances are more evident within fine-scale units. Analyses based on a single administrative unit or fixed grid can smooth local degradation patches, while the estimated relative importance of explanatory variables can change because of the modifiable areal unit problem (Liu et al., 2023b). Recent studies suggest that natural factors often retain strong explanatory power at local scales, whereas socioeconomic factors become more prominent at coarser scales (Xu et al., 2025a). Moreover, ecological quality generally has nonlinear relationships with natural conditions, human activities, and locational characteristics (Xia et al., 2026). An interpretable machine-learning framework combining Extreme Gradient Boosting (XGBoost) and SHapley Additive exPlanations (SHAP) can estimate variable contributions to model predictions and characterize nonlinear response patterns and interaction effects (Xu et al., 2025b; Zhao et al., 2026). These considerations motivate the integration of multiscale spatial analysis with interpretable machine learning to examine how ecological patterns and modeled associations vary with analytical scale.

The Hangzhou Bay Urban Agglomeration (HBUA), located on the eastern coast of China, is an important component of the Yangtze River Delta urban agglomeration. It is characterized by a high level of urbanization, intensive industrial activity, and strong land–sea interactions. The region encompasses diverse geographic units, including river-network plains, estuaries and bays, hills and mountains, coastal wetlands, and densely built-up areas, resulting in marked differences in baseline ecological conditions. Urban expansion, transport infrastructure development, industrial concentration, and ecological protection and restoration have long coexisted in the region (Sun et al., 2025; Xu et al., 2024). The HBUA plays an important role in the integrated development of the Yangtze River Delta, the concentration of advanced manufacturing, and China’s coastal opening-up. However, rapid urbanization and intensive land development also pose challenges to ecosystem stability and ecological connectivity (Huang et al., 2024). The combination of complex physical gradients and intensive human disturbance makes the HBUA a suitable region for evaluating long-term RSEI monitoring, multiscale ecological-quality assessment, and nonlinear associations with potential explanatory variables (Sun et al., 2022).

This study examines how the choice of spatial analysis unit changes conclusions drawn from a common ecological-quality dataset in the HBUA. A 30 m annual RSEI series for 2000–2024 is summarized using 1, 3, 5, and 10 km grids and county-level units. Separate XGBoost–SHAP models use pooled observations from six benchmark years (2000, 2005, 2010, 2015, 2020, and 2024). The study asks whether spatial aggregation changes (1) the mapped proportions of ecological improvement, stability, and degradation; (2) the ranking and relative predictive contributions of explanatory variables; and (3) the nonlinear associations and model interactions represented at each support. Its main contribution is the explicit comparison of ecological assessment and model interpretation across spatial supports. The comparison distinguishes conclusions that remain consistent across units from those that depend on aggregation, providing a basis for linking regional monitoring with local verification.

Study area

HBUA is located on the southern flank of the Yangtze River Delta, China, extending approximately from 118°20′ to 123°20′ E and from 28°50′ to 32°10′ N. The study area comprises Shanghai Municipality, Suzhou in Jiangsu Province, and Hangzhou, Ningbo, Huzhou, Jiaxing, Shaoxing, and Zhoushan in Zhejiang Province (Figure 1). It is an important component of the Yangtze River Delta urban agglomeration and represents a typical bay-area urban agglomeration characterized by a high level of urbanization, intensive industrial activity, and strong land–sea interactions.

The study area has diverse physical conditions and encompasses river-network plains, estuaries and bays, hills and mountains, coastal wetlands, and islands. Elevation generally decreases from the southwest to the northeast. The region has a subtropical monsoon climate with relatively abundant moisture and heat, and precipitation is concentrated mainly in spring and summer. Complex topographic gradients, extensive river networks, and coastal environments jointly shape the spatial heterogeneity of ecological quality.

Over the past two decades, the HBUA has experienced rapid expansion of construction land, densification of transport networks, spatial reorganization of industrial activities, and population concentration. Intensive human activities have altered land-use patterns and land-cover conditions, affecting vegetation cover, land surface temperature, surface moisture and dryness, and ecological connectivity. Concurrently, ecological restoration, wetland conservation, and urban green-space development have contributed to improvements in ecological quality. The region therefore provides a representative case for long-term RSEI monitoring, multiscale identification of ecological-quality patterns, and analysis of nonlinear associations with potential drivers.

Data sources and processing

Data sources

This study integrated remote sensing, geographic environmental, and socioeconomic spatial data to establish a database for ecological-quality assessment and explanatory-variable analysis. All data used to calculate RSEI were obtained through GEE, including Landsat 5/7/8/9 Collection 2 Level-2 surface reflectance and land surface temperature products and land use data from the China Land Cover Dataset (CLCD) for 2000–2024. A digital elevation model (DEM), precipitation, potential evapotranspiration, nighttime light intensity, and road-network data were also included. Table 1 lists the data sources and spatial resolutions.

Multiscale spatial aggregation

The five spatial supports were selected to span a practical range of aggregation and compare regular grids with administrative reporting units. Square grids with side lengths of 1, 3, 5, and 10 km have nominal uncut areas of 1, 9, 25, and 100 km², respectively (Figure 2). The 1 km grid retains relatively fine spatial contrasts in the 30 m RSEI data; the 3 and 5 km grids provide intermediate aggregation levels; and the 10 km grid summarizes broader spatial contrasts. County units provide a governance-relevant comparison with irregular boundaries and unequal areas. No scale was assumed to be optimal or universally applicable. The grid comparison mainly examines changes in support size, whereas comparison with counties combines differences in size, shape, zoning, and sample size (Jelinski and Wu, 1996).

The regular grids were generated in ArcGIS and intersected with the HBUA terrestrial boundary; county divisions supplied the administrative units (Table 2). RSEI was calculated at the pixel level before spatial aggregation; PCA was not refitted within individual grids or counties. Each unit summarizes the valid pixels in the same annual RSEI image.

Processing of explanatory variables

The potential explanatory variables represented topography, climate, urbanization pressure, and locational conditions (Huang et al., 2024; Xia et al., 2026). The topographic factors included elevation (DEM) and aspect (ASP), which characterized the terrain background. The climatic factors comprised precipitation (PRE) and potential evapotranspiration (PET), representing regional moisture and thermal conditions. Land use intensity (LUI) and night-time light intensity (NLI) were used to characterize urbanization pressure, including development activity and the intensity of human activity. Distance to road (DRD) represented transport accessibility and local disturbance.

The explanatory variables were summarized for 2000, 2005, 2010, 2015, 2020, and 2024 and matched to RSEI by year and spatial unit. PRE, PET, LUI, NLI, and DRD used the corresponding benchmark-year source layers; DEM and ASP were treated as static terrain attributes. Continuous variables were summarized using zonal means. Aspect angles were converted to sine and cosine, averaged within each unit, and transformed back into a circular mean angle (Hodgson and Gaile, 1996). This avoids arithmetic averaging across 0°/360°. LUI was computed from CLCD using weights of 7 for cropland, 10 for impervious surfaces, 2 for forest, shrubland, grassland and wetland, and 1 for snow/ice and bare land (Brown and Vivas, 2005). DRD is represented by the mean of pixel-level road-distance values within each statistical unit.

Before model development, variance inflation factor (VIF) tests were performed on the candidate variables to reduce the potential influence of multicollinearity on model fitting and interpretation (O Brien, 2007). The final XGBoost–SHAP models included seven factors: DEM, ASP, PRE, PET, LUI, NLI, and DRD (Table 3; Figure 3). The VIF values of all factors were below 5 at every analytical scale, indicating no problematic multicollinearity under the adopted criterion.

Methods

RSEI model

The RSEI integrates four remote sensing components—greenness, wetness, dryness, and heat—to characterize surface ecological conditions. Principal component analysis (PCA) was employed to reduce the influence of subjective weighting on the composite assessment (Hu and Xu, 2018; Wang et al., 2025; Xu, 2013). In this study, annual RSEI data for the HBUA from 2000 to 2024 were generated on the GEE platform using imagery acquired by Landsat 5 TM, Landsat 7 ETM+, Landsat 8 OLI/TIRS, and Landsat 9 OLI-2/TIRS-2.

First, valid Landsat observations acquired from June to September of each year were selected (Table 4). Clouds, cloud shadows, snow, and radiometrically saturated pixels were masked using the QA_PIXEL and QA_RADSAT bands. Open-water pixels were then excluded using the annual China Land Cover Dataset (CLCD). The remaining observations were composited using the annual median, after which NDVI, WET, NDBSI, and LST were calculated according to the formulas presented in Table 5.

In Table 4, ρ B l u e, ρ G r e e n, ρ R e d, ρ N I R, ρ S W I R 1 and ρ S W I R 2 denote the scaled Collection 2 Level-2 surface reflectance in the blue, green, red, near-infrared, shortwave-infrared 1, and shortwave-infrared 2 bands, respectively. D N S T _ B 6 and D N S T _ B 10 denote the quantized pixel values of the Level-2 surface-temperature bands for Landsat 5/7 and Landsat 8/9, respectively. LST is expressed in degrees Celsius (°C).

Second, each indicator was clipped and normalized to [0, 1] using common bounds over 2000–2024: the minimum annual 1st percentile and maximum annual 99th percentile. The same bounds were applied to every year.

Following the common-basis rationale of previous studies (He et al., 2025b; Xue et al., 2024), a fixed reference PCA was adopted. The normalized annual images were averaged using available valid observations to create a multiyear reference image. PCA was fitted once to the covariance matrix of this spatially mean-centered image (Jolliffe and Cadima, 2016). Annual PC1 scores were computed by centering the normalized annual indicators with the same reference means and applying the fixed, ecologically oriented loadings: +0.664386 (NDVI), +0.173277 (WET), −0.721935 (NDBSI), and −0.085886 (LST). PC1 explained 83.0135% of the reference-image variance.

Finally, Annual PC1 were calculated using common bounds defined by the minimum annual 1st percentile and maximum annual 99th percentile. Higher final RSEI values indicate more favorable composite surface conditions. Equal-interval classes were poor (0–0.2), fair (0.2–0.4), moderate (0.4–0.6), good (0.6–0.8), and excellent (0.8–1.0) (Hasan et al., 2025).

Cross-sensor agreement was evaluated over the complete study-area images using independent June–September composites for Landsat 5/7 (2000, 2005, 2010, and 2011), Landsat 7/8 (2013, 2015, 2020, and 2021), and Landsat 8/9 (2022–2024). The complete-image comparisons included 59,353,351–68,910,848 common-valid pixels per year (85.75–99.56% of the fixed terrestrial grid). For Landsat 5/7, r was 0.820–0.911, and RMSE was 0.105–0.157. For Landsat 7/8, r was 0.910–0.940 and RMSE was 0.113–0.138. For Landsat 8/9, r was 0.933–0.971 and RMSE was 0.070–0.100 (Table 6). Overall, RSEI estimates from different sensors show a degree of comparability in characterizing relative spatial patterns of ecological quality within the study area.

Hodrick–Prescott filter

The Hodrick–Prescott (HP) filter (Hodrick and Prescott, 1997) was applied to the regional mean RSEI series for 2000–2024 to distinguish its long-term trajectory and stage-specific characteristics from short-term anomalous fluctuations. The observed series was decomposed into a long-term trend component and a short-term fluctuation component (Equation 1):

where R S E I t is the regional mean RSEI in year t, τ t is the long-term trend estimated by the HP filter, and c t is the short-term fluctuation component. The smoothed trend was obtained by minimizing the following objective function (Equation 2):

The first term measures the deviation of the trend component from the original RSEI series, whereas the second controls the smoothness of the trend. The smoothing parameter is denoted by λ. A larger λ produces a smoother trend, whereas a smaller value allows the trend to follow the original RSEI series more closely. Because annual observations were used, λ was set to 6.25 according to the frequency-adjustment method (Ravn and Uhlig, 2002).

Theil–Sen slope estimation and Mann–Kendall test

The Theil–Sen median slope estimator (Sen, 1968) was used to calculate the rate of change in pixel-level RSEI. For any pixel p, the trend slope was defined as the median of the slopes calculated from all pairs of observation years (Equation 3):

where β is the annual rate of change in RSEI. A value of β >0 indicates an increasing trend, β <0 indicates a decreasing trend, and β =0 indicates no trend. To prevent minor fluctuations from being interpreted as directional change, values satisfying - 0.0005 ≤ β < 0.0005 were classified as stable. The Z statistic from the nonparametric Mann–Kendall test (Mann, 1945) was subsequently used to evaluate trend significance. Pixel-level trends were divided into five categories according to the direction of the Theil–Sen slope and the Mann–Kendall test result (Peng, 2024) (Table 7).

XGBoost–SHAP analysis

XGBoost (Chen and Guestrin, 2016) and SHAP (Lundberg et al., 2020; Lundberg and Lee, 2017) were used to examine associations between RSEI and the selected explanatory variables across spatial supports. The response was the RSEI level for a spatial unit in a benchmark year. Each row represented one spatial unit in 2000, 2005, 2010, 2015, 2020, or 2024. The pooled models therefore combine between-unit and between-year variation without identifying separate within-unit temporal effects. Rows with nonfinite values or the −9999 missing-value code in the response or any selected predictor were excluded. The resulting sample sizes were 363,395, 43,819, 16,707, 4,697, and 426 unit-year records at the 1 km, 3 km, 5 km, 10 km, and county supports. XGBoost generated predictions by integrating multiple regression trees (Equation 4):

where y ^ i is the predicted RSEI for sample i, X i represents its explanatory variables, f k is the (k)th regression tree, and K is the number of trees. At each spatial scale, the data were divided into training, validation, and test sets at proportions of 64%, 16%, and 20%, respectively, using a random seed of 42. The model employed the squared-error regression objective, with root mean squared error (RMSE) used as the metric for monitoring model training. The maximum number of boosting rounds was set to 900, the maximum tree depth to 4, and the learning rate to 0.05. Training was terminated early if the validation RMSE did not improve for 60 consecutive rounds, and predictions used all trees retained when early stopping terminated training. Additional settings were subsample = 0.8, colsample_bytree = 0.8, L2 regularization = 5, L1 regularization = 1, and minimum child weight = 3.

SHAP was used to quantify the contribution of each explanatory factor to the model prediction (Equation 5):

where φ 0 is the model baseline value, φ ij is the SHAP value of factor j for sample i, and M is the number of explanatory factors. A positive SHAP value indicates that the corresponding factor increases the predicted RSEI relative to the model baseline for the given observation, whereas a negative value indicates that it decreases the prediction. The relative contribution of factor j was calculated from its mean absolute SHAP value (Equation 6):

where C j is the relative contribution of factor j, and n is the number of samples. SHAP dependence plots, partial dependence plots (PDPs), and SHAP interaction values were then used to characterize nonlinear responses, identify thresholds, examine model interactions, and compare these patterns across spatial scales (Wu et al., 2025).

For each predictor, the PDP was evaluated at 45 equally spaced values between its 1st and 99th percentiles using at most 600 randomly sampled background rows. PDP means and the bounds from 120 background-row resamples were smoothed by a centered five-point moving average. The reported reference position maximizes the absolute numerical PDP gradient after excluding the outer 10% of evaluation-grid positions and at least two positions at each edge.

RSEI-Based ecological quality assessment

Figure 4 presents the spatial distribution of RSEI in six representative years. The regional mean RSEI values in 2000, 2005, 2010, 2015, 2020, and 2024 were 0.5959, 0.5280, 0.5654, 0.6091, 0.6148, and 0.6382, respectively, indicating an initial decline followed by recovery. RSEI decreased by 11.38% from 2000 to 2005 and subsequently recovered. By 2024, it was 0.0423 (7.10%) higher than in 2000. Across the complete annual series, the mean RSEI was 0.5854 and ranged from 0.5280 to 0.6471. The minimum occurred in 2005 and the maximum in 2023. Although RSEI decreased slightly from 2023 to 2024, the 2024 value remained relatively high within the study period.

The four RSEI components exhibited distinct interannual dynamics (Figure 5a). LST fluctuated between 0.744 and 0.849 and reached 0.833 in 2024, exceeding its 2000 value of 0.764. NDBSI reached its maximum of 0.512 in 2005 and then declined overall to 0.385 in 2024. NDVI decreased to 0.693 in 2008 and subsequently recovered with fluctuations, reaching 0.758 in 2024. WET decreased to 0.315 in 2010 and generally increased thereafter, reaching 0.409 in 2023 and 0.403 in 2024. Given the PC1 loading directions, the later recovery in RSEI coincided with increases in NDVI and WET and a decrease in NDBSI, whereas LST remained negatively oriented in the index. The HP filter further identified a nonmonotonic long-term trajectory (Figure 5b). The trend component initially declined and then increased, with an overall increase of 8.85% during the study period. The short-term fluctuation component ranged from −0.0269 to 0.0166, with the largest negative deviation occurring in 2012 and the largest positive deviation in 2017.

Spatiotemporal changes in RSEI

Spatial distribution

Across the five spatial scales, RSEI exhibited a persistent gradient characterized by high values in the southwest and relatively low values in the northeast and densely urbanized areas along Hangzhou Bay (Figure 6; Supplementary Figure S1). Higher values occurred mainly in the southwestern hills and mountains, contiguous forested areas, and locations with relatively continuous ecological spaces. Lower values were concentrated along the high-intensity development corridor connecting Shanghai, Suzhou, Jiaxing, Hangzhou, and Ningbo, as well as along urban fringes and in some areas bordering the bay.

The overall spatial gradient persisted across scales, but local heterogeneity progressively diminished with increasing spatial aggregation. The 1 and 3 km grids retained the interspersed high- and low-value patches associated with urban fringes, transport corridors, islands, and coastal areas. In contrast, the 10 km grid and county-level units emphasized regional mean conditions and contained fewer units with extremely high or low values.

Changes in RSEI

The full-period Theil–Sen trends and Mann–Kendall test results indicated that improvement predominated at all five spatial scales between 2000 and 2024 (Figure 7). On an area-weighted basis, the combined proportions classified as slight or significant increases were 75.34%, 76.91%, 77.71%, 79.49%, and 84.25% at the 1 km, 3 km, 5 km, 10 km, and county scales, respectively. The corresponding combined proportions of slight and significant decreases declined from 20.64% to 14.35%, while stable areas accounted for 1.40%–5.42%. Significant increases represented 18.00%–21.03% of the study area across scales and occurred mainly in the southwestern hilly forestlands, contiguous ecological spaces in the south, and some coastal and island areas. Slight and significant decreases were more prevalent in the northern plain, densely urbanized areas along the bay, and major urban-expansion corridors. As the analytical unit increased from the 1 km grid to the county scale, the proportion of slight increases rose from 57.35% to 63.44%, whereas the proportion of significant decreases declined from 2.11% to 0%. Thus, the long-term predominance of improvement was consistent across scales, although spatial aggregation reduced the mapped extent of local degradation and extreme change.

Stage-specific RSEI changes were also scale dependent (Figure 8). During 2000–2005, decreases predominated at every scale, and the area with a decrease of more than 2%—classified as degradation—ranged from 58.36% to 70.32%, with the highest proportion at the county scale. The direction of change reversed during 2005–2010 and 2010–2015. The area with increases of more than 2%—classified as improvement—accounted for 59.72%–69.18% and 68.41%–87.71%, respectively; the highest value in the latter interval occurring at the county scale. The regional pattern therefore shifted from widespread decline in the early period to extensive recovery in the middle periods, although the estimated area of improvement varied among spatial units.

The strongest scale differentiation occurred during 2015–2020. At the 1 km scale, degradation, stability, and improvement accounted for 33.52%, 27.22%, and 39.26%, respectively; stability was defined as a change between −2% and 2%. At the county scale, the corresponding proportions were 22.52%, 60.85%, and 16.64%. Relative to the 1 km results, county-level aggregation increased the stable proportion by 33.63 percentage points and reduced the improved proportion by 22.63 percentage points. Spatial aggregation therefore altered the dominant change category: fine-scale grids retained a mosaic of degradation and improvement patches, whereas coarser grids and county-level units more often classified opposing local changes as stable.

Multiscale association analysis

Scale-dependent factor contributions

The XGBoost models showed strong predictive performance at all five spatial scales (Figure 9). The test-set R² values at the 1 km, 3 km, 5 km, 10 km, and county scales were 0.830, 0.867, 0.858, 0.843, and 0.931, respectively, while the corresponding RMSE values were 0.083, 0.067, 0.066, 0.067, and 0.051. The county-scale model had the highest test-set R² and the lowest RMSE. Across the grid scales, test-set R² did not vary monotonically from 3 to 10 km.

The mean absolute SHAP values identified LUI as the most important factor associated with RSEI across scales, with a mean relative contribution of 51.46% across the five scales (Figure 10). Its relative contributions at the 1 km, 3 km, 5 km, 10 km, and county scales were 55.34%, 53.52%, 55.43%, 48.43%, and 44.59%, respectively, ranking first among all factors at every scale. PRE, NLI, DEM, and DRD were secondary factors, with mean relative contributions across the five scales of 14.17%, 9.89%, 9.86%, and 8.14%, respectively.

The other variables also exhibited scale-dependent contributions. PRE increased from 12.44% at the 1 km scale to 15.67% at the county scale, and NLI increased from 6.80% to 20.33%. DRD contributed 7.52%–10.00% across the 1–10 km grid scales but only 5.71% at the county scale. ASP and PET made comparatively small contributions.

The dependence plots revealed nonlinear modeled associations across scales (Supplementary Figure S2). Higher LUI was associated with lower predicted RSEI at every support. Figure 11 summarizes the descriptive thresholds for LUI, PRE, NLI, and DRD across the five spatial supports. The thresholds were 7.54–7.71 for LUI across the grid supports and approximately 8.25 for counties. Corresponding thresholds were approximately 1208–1266 mm and 1305 mm for PRE, and approximately 17.25–19.31 km and 5.59 km for DRD. For NLI, they were 3.92, 3.67, 3.55, 3.34, and 5.02 at the five supports, respectively.

Interaction effects

The SHAP interaction ranking changed with spatial scales (Figure 12). DEM–LUI ranked first at 1 km, 3 km, 5 km, and 10 km, with mean absolute interaction values of 0.0320, 0.0307, 0.0138, and 0.0113, respectively. At the county scale, LUI–NLI ranked first at 0.0034, PRE–PET, PRE–LUI, DEM–PET, DEM–NLI, and PET–DRD also had relatively high magnitudes at some scales. Thus, the highest-ranked interaction shifted from a topography–land use pair at the grid scales to a land use–urbanization pair at the county scale.

Long-term trajectory and spatial differentiation of ecological quality

Ecological quality in the HBUA did not improve monotonically between 2000 and 2024 but followed a stage-specific trajectory characterized by an early decline, a mid-period transition, and sustained recovery during the later period. The HP-filtered results further showed that this long-term recovery trend was superimposed by pronounced interannual fluctuations. Ecological quality improvement over the study period should therefore be interpreted as an overall long-term tendency rather than continuous and synchronous improvement across all locations and years.

Changes in the four RSEI components were consistent with this nonmonotonic trajectory. During the later study period, NDVI and WET generally increased, whereas NDBSI declined, indicating vegetation recovery, greater surface wetness, and a lower composite dryness signal associated with bare and impervious surfaces. LST remained relatively high, however, suggesting that the recovery of greenness and wetness did not fully offset thermal pressure under rapid urbanization. These concurrent patterns imply that improved composite ecological quality in a densely urbanized agglomeration can coexist with persistent thermal pressure. A vegetation indicator alone is therefore insufficient to represent the overall surface ecological condition. Previous remote sensing studies likewise indicate that multidimensional surface indicators provide a more complete representation of urban ecological conditions than greenness alone and that urban green-space cooling depends on size, quality, spatial configuration, and connectivity (Li et al., 2024; Yu and Fang, 2023).

Spatially, the HBUA maintained a relatively stable RSEI pattern characterized by high values in the southwest and low values in the northeast and densely urbanized areas along Hangzhou Bay. The southwestern hills and mountains have extensive forest cover and relatively continuous ecological spaces. Their topographic conditions may constrain intensive construction and help maintain vegetation growth and surface hydrothermal regulation. By contrast, the intensively developed corridor connecting Shanghai, Suzhou, Jiaxing, Hangzhou, and Ningbo contains extensive impervious surfaces, transport infrastructure, and concentrated population and industrial activities. Green-space fragmentation and heat accumulation may jointly generate compound ecological pressure in this corridor (Zhang et al., 2025). Areas along the bay and around the estuary are concurrently affected by urban expansion, wetland changes, and land–sea interactions, resulting in greater spatial complexity in ecological-quality changes.

RSEI represents a composite surface ecological condition characterized by greenness, wetness, dryness, and heat; it should not be considered equivalent to regional ecosystem integrity. A high RSEI does not necessarily indicate similarly high levels of biodiversity, habitat connectivity, air quality, or nearshore water quality. For example, some areas with extensive vegetation cover may still experience landscape fragmentation, simplified community structure, or disruption of ecological corridors. Previous research has also emphasized that the ecological benefits of green infrastructure depend on connectivity, multifunctionality, and adaptive management rather than increased green cover alone (Wang et al., 2024). Accordingly, the spatial differences identified in this study are more appropriately used to delineate zones for protection, restoration, and priority diagnosis.

Scale dependence of RSEI change detection

One of the principal findings of this study is that spatial aggregation not only reduced local variation in the RSEI distribution but also systematically altered the classification of ecological degradation, stability, and improvement. Although ecological improvement predominated at all five scales over 2000–2024, the area proportions assigned to individual change categories differed among scales. Scale differentiation was greatest during 2015–2020: relative to the 1 km scale, the county scale increased the stable proportion by 33.63 percentage points and reduced the improved proportion by 22.63 percentage points (Figures 13a–c).

This result indicates that spatial aggregation does more than remove random noise. When improvement and degradation patches coexist within the same coarse-scale unit, changes in opposing directions may cancel each other in the mean and consequently be classified as stable. Conversely, when one direction of change clearly predominates within an aggregated unit, spatial aggregation may reinforce that change category (Oshan et al., 2022). A coarser scale therefore does not necessarily increase the stable proportion. The direction of its influence depends on the composition, spatial configuration, and relative proportions of local change patches within each aggregated unit. This dependence of assessment results on spatial-unit size and zoning scheme is a typical manifestation of the MAUP (Ju et al., 2021).

Differences between regular grids and county-level administrative units further indicate that the scale effects identified here include changes in both spatial-unit size and zoning scheme. Regular grids have relatively consistent shapes and areas and can represent continuous spatial gradients more directly. County-level units, however, vary considerably in area, shape, land-use composition, and internal heterogeneity. Their boundaries are primarily designed for administrative management and do not necessarily correspond to the boundaries of ecological processes. Consequently, county-level means can obscure local degradation along urban expansion fronts, transport corridors, rivers, and coastal development zones, while accentuating regional differences in overall ecological conditions and urbanization levels.

Cross-scale changes in factor contributions further demonstrate the sensitivity of model interpretation to spatial support (Figure 13d). Relative to the 1 km scale, differences in factor contributions were limited at the 3 km scale but increased at 5 and 10 km. The largest restructuring occurred at the county scale, where the relative contribution of NLI increased, while those of LUI, DEM, PET, and DRD decreased. This pattern is consistent with the smoothing of fine-scale land use and road-distance variation through aggregation and the greater contrast in county-level urbanization and hydroclimatic conditions.

The results therefore do not support selecting a single “optimal scale” independently of the research objective. The 1–3 km grids preserve variation associated with urban fringes, transport corridors, and local ecological patches and can support hotspots screening. The 5–10 km scales suppress fine-scale variation and can support assessment of landscape-scale gradients and cross-boundary patterns, but they do not establish the true operating scales of ecological processes. The county scale corresponds to common policy and territorial-governance units and is useful for regional comparison and priority screening. The value of multiscale analysis lies in distinguishing the observation scale, the hypothesized process scale, and the policy-implementation scale (Saberski et al., 2024).

Scale-dependent associations and nonlinear interactions

The fitted models explained substantial RSEI variation under the reported validation design. The county model’s higher random-split R² does not establish that counties represent ecological processes more accurately: aggregation changes the target variance, sample composition, and predictor relationships (Jiang et al., 2024), while random splitting retains spatial and repeated-unit dependence between subsets.

LUI had the highest relative contribution at every scale. This finding is consistent with previous research (Zhao et al., 2026). Within the available predictor set, land-use composition therefore provided information used extensively by the fitted models. This does not establish that land use has a greater causal influence than precipitation: predictors differ in native resolution and spatial variability, and land cover is closely related to the greenness, dryness, and heat components of RSEI. The share assigned to NLI rose from 6.80% at 1 km to 20.33% at the county support, whereas DRD had a smaller share for counties. Such differences are consistent with aggregation altering the spatial contrasts represented by the models (Mai et al., 2025).

The nonlinear dependencies and interaction results further indicate that the factor contributions were not fixed or purely additive. The high-gradient intervals of LUI, NLI, PRE, and DRD shifted among scales, reflecting changes in model-sensitive ranges under different spatial supports. Threshold analyses should therefore explicitly specify the corresponding spatial scale rather than treating the identified thresholds as universally applicable (Saberski et al., 2024). In this study, DEM–LUI was the dominant interaction pair across the 1–10 km grid scales, suggesting that topographic conditions may modify the association between land-use intensity and ecological quality. At the county scale, the leading pair shifted to LUI–NLI, indicating that land use intensity and broader urbanization activity were jointly represented in the county-level model. Similar nonlinear associations and interaction patterns have been reported for the Yangtze River Economic Belt and other urban agglomerations (Liu et al., 2023a; Wei and Huang, 2025).

Implications for multiscale ecological governance

The findings can inform a scale-nested ecological monitoring framework that links regional screening, cross-boundary pattern assessment, hotspot diagnosis, and site-level verification. This framework aligns remote sensing outputs with different governance levels while retaining the need for local validation before intervention.

The county scale can be used to compare overall ecological conditions, urbanization pressure, and long-term change directions and to screen administrative units for further assessment or cross-regional coordination. The 5–10 km scales can help identify hill–plain transition zones, river-network and wetland systems, ecological corridors along Hangzhou Bay, and candidate cross-boundary patterns. The 1–3 km scales can further locate urban expansion fronts, road-adjacent areas, concentrated development zones, and local degradation patches. Planning parcels, high-resolution imagery, and field observations should then be integrated with the 30 m RSEI results before specific intervention boundaries are delineated. These recommendations are consistent with multiscale ecological-security research in the Yangtze River Economic Belt (Wei and Huang, 2025).

The observed spatial patterns suggest several hypotheses for differentiated governance, subject to verification with local planning and field data. In the southwestern hilly forestlands and contiguous ecological spaces, assessment may prioritize habitat continuity, fragmentation risk, and construction expansion. In the densely urbanized corridor connecting Shanghai, Suzhou, Jiaxing, Hangzhou, and Ningbo and along urban expansion fronts, follow-up assessment may examine impervious-surface control, blue–green networks, thermal exposure, and ecological connectivity. Wetlands along the bay, rivers crossing county boundaries, and mountain corridors also warrant cross-jurisdictional monitoring at landscape scales rather than isolated interpretation of county means.

Ecological-governance priorities should not be determined from RSEI alone. Where independent exposure and service data confirm high population density, high heat exposure, insufficient green-space provision, and low RSEI, such areas may warrant greater attention than areas with similar RSEI but lower population exposure. Previous research has shown that the cooling benefits of urban green spaces are unevenly distributed across populations and locations (Li et al., 2024). Future governance assessments should therefore integrate RSEI with population distribution, heat exposure, green-space accessibility, habitat connectivity, and water-quality indicators. Such integration can reduce the risk that improvement in a greenness-related indicator obscures other ecological costs (Wang et al., 2024).

Limitations and future research prospects

This study has several limitations. First, RSEI represents greenness, wetness, dryness, and heat rather than biodiversity, habitat connectivity, air quality, or ecosystem integrity. Second, explanatory models pool six benchmark years, rather than a complete annual covariate panel, and do not isolate within-unit temporal effects. Third, source-resolution differences, filled factor values, small clipped units, and no temporal holdout constrain generalization. No spatial block or spatial-group validation or temporal holdout was performed. Random splitting may place spatially dependent observations and repeated observations of the same unit in different subsets, potentially inflating predictive performance. Fourth, the analysis focused primarily on changes in spatial-unit size and did not systematically evaluate the zoning effects caused by shifts in grid origin. Finally, the analysis compared the relative importance of the factors but did not conduct causal inference.

Future work could integrate Sentinel and high-resolution thermal observations, biodiversity and connectivity measures, and independent environmental monitoring. Additional spatial and temporal holdouts, predictor-resolution sensitivity checks, and panel designs separating within-unit change from between-unit differences would improve interpretation. Spatial panel models, natural experiments, and causal-inference designs would be needed to evaluate the effects of specific interventions.

This study constructed an annual RSEI dataset for the HBUA from 2000 to 2024 to compare ecological-quality changes across scales and modeled associations using six benchmark-year datesets. The principal conclusions are as follows:

Ecological quality in the HBUA exhibited a nonmonotonic long-term trajectory. The regional mean RSEI declined substantially from 2000 to 2005 and subsequently recovered overall, reaching 0.6382 in 2024. This was 0.0423 (7.10%) higher than the 2000 value. Spatially, RSEI maintained a relatively stable pattern of high values in the southwest and low values in the northeast and densely urbanized areas along Hangzhou Bay. The later improvement coincided with increases in NDVI and WET and a decrease in NDBSI, while persistently high LST indicated that thermal pressure did not decline concurrently.

Ecological quality changes exhibited pronounced scale dependence. As the analytical unit increased from the 1 km grid to the county scale, the proportion of slight increases rose from 57.35% to 63.44%, while the proportion of significant decreases declined from 2.11% to 0%. Administrative-unit averages may classify coexisting improvement and degradation patches that offset one another within the same unit as stable.

Model-derived importance of the associated factors changed with spatial support. LUI had the largest predictive contribution at every scale. PRE, DEM, NLI and DRD were secondary factors. The contribution of NLI reached 20.33% at the county scale, whereas the share assigned to DRD decreased. Thus, regional contrasts in urbanization were represented more strongly at the county scale, while road-distance associations were represented more strongly at local and intermediate scales.

Nonlinear associations and interactions were also scale specific. Fitted response shapes and thresholds varied with aggregation. DEM–LUI was the leading interaction in the grid models, whereas LUI–NLI ranked first for counties.