The spatiotemporal analysis model includes three parts: (1) Thei-Sen (Sen) is a method used for trend analysis (Feng et al., 2020), enabling the analysis and description of the trends in ESs from 2000 to 2022; (2) Mann–Kendall (MK) is a significance testing approach (Feng et al., 2020; Feng X. et al., 2023), and it can eliminate a few outliers with a good adaptability (Hensel and Frans, 2006); and (3) Hurst is a method for portraying the information dependence (Zhang and Jin, 2021), which helps determine whether the evolution of ESs follows continuous patterns through via time series. It is in Hurst that we use the rescaled range analysis method (R/S) (Jiang et al., 2017) to uncover the evolutionary features. In this study, we investigated the spatiotemporal evolution of ESs by coupling Sen, MK and Hurst. The rationale is shown in Figure 3.

Where y _i and y _j stand for the ESs of monitoring years i and j, respectively, and θ represents the trends in ESs. D and Z indicate the test statistics and the standard test statistics, respectively; n stands for the samples; it is y _i −y _j > 0 that Sign = 1, and vice versa, Sign = −1. In Hurst, X(t, τ) denotes the sequences of cumulative deviations; R(τ) and S(τ) stand for the range sequences and the standard deviation sequences, respectively; Hurst is between 0 and 1; and it is H > 0.5 that ESs have a positive continuous trend with a persistence feature and vice versa. Moreover, we use spatial principal component analysis (SPCA) to address multiple variables, replacing original variables with a composite variable (Feng L. et al., 2023; Junttila and Laine, 2017). SPCA eliminates the need for artificially determined weights and lowers the subjectivity the evaluation of ESs. In this study, Pc _i denotes the principal information sequences i; and w _i indicates the weight of Pc _i . When the cumulative variance contribution is greater than 80%, it can represent most of the information.

Random forest (RF) is, based on statistical theory, a machine learning method (Breiman, 2001), which is by means of an average combination of multiple decision trees to obtain the final regression result (Marco et al., 2022). RF avoids overfitting with a good robustness. With the help of the mean decrease Gini (MDG), we can categorize the importance of indicators. Equation 1 is as follows: G I n = 1 − ∑ k = 1 K P n k 2 (1) where GI _n is the reduction of node impurity by the factors n; P _nk denotes the percentage of features in nodes n; and K is the number of nodes. In addition, with the help of the partial dependence function for single-factor regression, we describe the effect of dominant factors by a partial dependence plot (PDP). Equation 2 is as follows: F ∧ X S = E X S F ∧ X S , X C = ∫ F ∧ X S , X C d P X C (2) where Xs represents the dominant influencing factor and Xc stands for all other influencing factors. Moreover, the partial function is estimated by calculating the mean value of the training set with the aid of the Monte Carlo method.

Multiscale geographically weighted regression (MGWR), with a better spatial smoothing, is an enhancement of GWR, which can generate self-adaptive bandwidths at different scales (Yu et al., 2020). In addition, the weight of MGWR is expressed by a distance function between the observation point and the regression point, which aims to evaluate the importance to which the observations have an influence on the parameter estimates (Duan et al., 2021). Equation 3 is as follows: y i = ε i + ∑ j = 1 k β b w j u i , v i x i j (3) where y _i and x _ij indicates the explained variable and explanatory variable, respectively; i stands for the sample size; j is on behalf of the independent variable size; (u _i ,v _i ) and ε _i represents the spatial location and the error term, respectively; and bwj represents the bandwidth. Moreover, we use the corrected akaike information criterion (AICc) to evaluation the superiority in models; a lower AICc indicates a better model fit (Sisman and Aydinoglu, 2022); and the same applies to the bayesian information criterion (BIC) does. In the end, the variance inflation factor (VIF) is used to check for multicollinearity, and VIF<7.5 indicates weak collinearity among the variables.

As shown in Figure 4, ESs are at the middle level (Mean = 0.538), with strong temporal features and significant interactions. Figure 4A shows that ESs exhibit a fluctuating growth trend (k = 0.017, R ² = 0.175) from 2000 to 2022, which can be broadly categorized into four periods (T₁:2000–2005, T₂:2005–2012, T₃:2012–2016, T₄:2016–2022). ESs are better coordinated in the early stages (T₁, T₂ and T₃), undergoing an evolutionary process of decreasing, increasing, and then decreasing again. However, in the later stages (T₄), they gradually tend to become disordered poor consistency. As depicted in Figures 4B, C, ESs, with the first principal component contributing 82.52%, show a strong connection (Correlation = 0.724 and R ² = 0.538) with each other. However, En has a low correlation with other ESs, contributing less than 6.82%. As a result, En plays a minimal role in the actual evaluation of ESs. In summary, ESs show slight improvement, but their internal functions are at risk of imbalance. There are both competitive and cooperative relationships within ecosystems. If the disturbances from social systems are not considered, the steady-state transition in ecosystems is generally relatively slow. However, CZC has an inherent advantage in resource enrichment, which leads to a rapid aggregation of populations. This phenomenon can intensify competition between the systems, thereby accelerating the steady-state transition. Therefore, this internal dysfunction could be the result of high-intensity human activities.

The spatial evolution of ESs in Figure 5 is analyzed in conjunction with Figure 1. As shown in Figure 5A, spatially, Es has a clear north‒south discrepancy, characterized by a low level in Region N (including Region N₁, Region N₂ and Region N₃) and a high level in Region S (including Region S₁ and Region S₂). The low and middle-low levels of Es are primarily located in Region N, accounting for 36.5% of CZC; the high and middle-high levels of Es are mainly found in Region S, covering 42.3% of CZC; and the middle level of Es has a broken spatial distribution, occupying 21.2% of CZC. Firstly, Pr, Sc, Bi and Hq have a high similarity in their spatial distribution, with a significant north‒south difference. Specifically, the middle-low levels and below of Pr, Sc, Bi, and Hq (accounting for 40.9%, 54.9%, 35.2% and 59.2%, respectively) are positioned in Region N. Meanwhile, the middle levels and above of them (covering 59.1%, 45.1%, 64.8% and 40.8%, respectively) are found in Region S, exhibiting strong spatial heterogeneity, particularly in Region S₂. The climate zone of CZC transitions from temperate to tropical from north to south, resulting in the formation of multiple ecosystems (e.g., grassland ecosystems, cropland ecosystems, and woodland ecosystems) from high latitudes to low latitudes. Consequently, this cross-scale climate could be the reason for their north-south differentiation. Secondly, except for Region N₁, the spatial distributions of Cs and Ws are similar to those of Pr, Sc, Bi, and Hq. Specifically, Cs is lower in Region N₂ and Region N₃, but higher in other areas, while Ws (Mean = 0.771) is higher in all regions except for Region N₁. On the one hand, the spatial distributions of ESs rely on large-scale natural conditions (such as topography and climate). On the one hand, their spatial distributions are also affected by local environments and human activities. Hence, it is thought that the spatial differences between them could result from a combination of microclimates and human activities. Finally, En is relatively high (Mean = 0.799), with its high level covering 67.6% of CZC, but it lacks a clear spatial aggregation feature. It is known that CZC has a high landscape aesthetic value, which leads to a higher level of En.

As shown in Figure 5B, Es shows a slight improvement (k = 0.017), accompanied by localized degradation (occupying 26.7%). Spatially, the changes in the bipolar regions (Region N₁ and Region N₂ and Region S₂) of Es are relatively smooth, while those in the central regions (Region N₃ and Region S₁) are relatively drastic. Firstly, both Pr (k = 0.011) and Ws (k = 0.018) show improvement with a similar spatial distribution. The improvement areas are concentrated in Region N (accounting for 54.9% and 40.9%, respectively). However, in Region S, although the trends remain relatively stable, localized degradation occurs. Secondly, Bi (k = 0.010), Hq (k = 0.028) and En (k = 0.031) have an increasing trend with a significant north‒south divergence. The improvement areas are mainly localized in Region S (occupying 26.7%, 56.3% and 57.7%, respectively), exhibiting a consistent spatial distribution; in Region N, there is strong heterogeneity with a broken spatial distribution, and localized degradation areas are evident (accounting for 33.9%, 31% and 21.2%, respectively). Thirdly, Sc (k = 0.007) shows a relatively stable trend. Spatially, Region N of Sc is more stable than Region S. Eventually, Cs has a decreasing trend, with its degraded areas concentrated in the midland (including Region N₂, Region N₂ and Region S₁). However, Cs maintains a steady trend in the bipolar regions. In summary, ESs have an increasing trend as a whole, with significant spatial zoning characteristics. On the one hand, the oceanic climate (hydrothermal condition) guarantees that CZC is characterized by a good natural attribute; in the face of external disturbances, ESs show a high resilience, which could be the main reason for the increase in ESs. On the other hand, CZC spans multiple climatic zones with different phenology features; these climatic zones divide CZC into several regions, forming multiple subsystems, and the interregional disparity becomes wider under both natural and social influences; therefore, this zoning phenomenon could be a result of the intensified spatial processes of ESs.

As shown in Figure 6, Es is expected to experience slight improvement in the future, with a distinct zoning pattern. Spatially, Es is at risk of degradation in localized areas (Region N₁, Region N₂ and Region S₁), but has an improving phenomenon in Region S₂ (occupying 42.3%). Firstly, Pr (mean = 0.431), Cs (mean = 0.608) and Ws (mean = 0.417) are at risk of degradation (accounting for 45.1%, 53.5% and 35.1%, respectively), especially in Region N, and they have a lower risk in Region S with strong heterogeneity. In Region N, the natural conditions are relatively poor, and the ecosystem is dominated by cropland ecosystems. Hence, the ecosystem shows a lower resilience to external disturbances, which could be the reason why they have a higher risk in Region N. However, in Region S, this circumstance is opposite. For localized areas with higher risk, this phenomenon could be the result of human activities. Secondly, Sc (mean = 0.512) and En (mean = 0.496), in Region S, have a relatively steady trend with lower risk. In Region N, the terrain is relatively flat, and the landscape is rather homogenous, which can provide certain prerequisites for Sc and En and make them change indistinctively in the short term. However, these prerequisites have an obvious gap in Region S. Hence, it is thought that these prerequisites could be reasons for the north-south risk divergence in Sc and En. Finally, Bi (Mean = 0.523) and Hq (Mean = 0.538) are expected to improve in the future (occupying 45.1% and 56.3%, respectively). Bi and Hq are susceptible to human activities and climatic conditions. With the implementation of ecological management projects in recent years, ecological conservation achievements in some areas have begun to bear fruit. However, localized areas still face some ecological risks due to the specificities of their urban development. Therefore, Bi and Hq with a lower risk could be the result of ecological preservation. Overall, the risk evaluation of ESs can provide some reference for stakeholders. It is by real-time monitoring that stakeholders have timely access to information, which plays a positive role in promoting the sustainable development of CZC.

As shown in Table 2, MGWR, with better spatial smoothing capability, slightly outperforms GWR. This study tests the colinear effect of factors by the variance inflation factor and acquires the main influencing factors (VIF<7.5) in ESs. In MGWR, both AICc and BIC are lower than in GWR, though the differences are not significant. In conclusion, MGWR makes use of fewer parameters to approach the true values with a better explanatory ability.

As portrayed in Figure 8, terrain and soil have a stronger influence on ESs in terms of natural factors, while climate exerts a weaker influence. Sand (Mean = −0.143) and Tp (Mean = −0.024) have a negative impact on ESs, but Dem (Mean = 0.128), Slope (Mean = 0.135) and Tem (Mean = 0.097) show a positive impact on ESs. Dem (0.119–0.152) and slope (0.135–0.145) have a higher level in Region S (spatially occupying more than 41.4% and 36.3%, respectively) and have a lower level in Region N. However, Tem (0.010–0.156) has a lower level in Region S (accounting for 53.4%) and a higher level (0.156–0.218) in Region N. Adequate hydrothermal conditions are important for the improvement of ESs. Region N is characterized by inadequate hydrothermal conditions and lower ESs. Hence, ESs are spatially more susceptible to Tem in Region N. Complex terrain enhances the natural features of localized areas and somewhat limits high-intensity human activities, which is more conducive to the development of ESs in a favorable direction. Region S is characterized by better natural features and higher ESs. As a result, the terrain spatially shows a more positive influence on ESs in Region S. Sand (Mean = −0.143) has a negative effect on ESs, exhibiting a more significant influence (−0.150∼-0.141) in Region S (accounting for 60.4%) compared to Region N (−0.141∼-0.132, occupying 39.6%). CZC, with a high sand and gravel content, is adjacent to the sea and susceptible to the effects of soil erosion. As a result, ESs are more susceptible to the effects of Sand, especially in the geologically complex Region S. The effect of Tp on ESs shows a tendency toward 0, but it has a weak negative effect on Region S₁. Frequent extreme weather is prone to impact the structure and function of ecology, which can put a damper on the improvement of ESs. Hence, spatially, Tp exerts a weak negative effect on ESs.

For social factors, Cult (Mean = −0.316) and Urban (Mean = −0.280) have the strongest impact on ESs, followed by Pop (Mean = −0.128) and Pm (Mean = −0.171), while LST (Mean = −0.124) and Port (Mean = 0.075) have the least impact. Social factors (excluding Port) have a negative effect on ESs. Cult and Urban show lower levels (−0.314∼-0.295 and −0.272∼--0.271) in Region N (accounting for approximately 37.9% and 31.1%, respectively), but higher levels (−0.345∼-0.314 and −0.286∼-0.280, respectively) in Region S (accounting for 38.0% and 51.7%, respectively). The agglomeration feature of Cult and Urban can reflect the stability of ecological structure. For a long time, Region N is dominated by the ecosystem of city and cropland with lower ESs, so it had a stable ecological structure. However, the opposite is true for Region S. Thus, ESs in Region S are more susceptible to perturbation. Pop and Pm exhibit opposite spatial distribution patterns. Pop has a higher level (−0.160∼-0.126) in Region S (occupying 51.8%) but a lower level (−0.126∼-0.068) in Region N (accounting for 48.2%), while the opposite is true for Pm (−0.214∼-0.188 and −0.188∼-0.140 accounting for 48.3% and 51.7%, respectively). Region N exhibits higher air pollution, and Region S has a higher population density. Therefore, spatially, ESs are vulnerable to Pm in Region N and to Pop in Region S. LST and Urban exhibit similar spatial distributions; however, LST has a slightly positive influence on ESs in Region N. Region N experiences poorer light and heat conditions compared to Region S, while LST creates favorable temperature conditions for crop growth to some extent. This could be the reason why LST has a positive effect in Region N. Moreover, Port is close to 0 and has a slight negative influence on ESs. Port can disturb the ecological environment, but due to their small size and low level, the disturbances are localized. As a result, spatially, the effect of Port on ESs is not significant.

In Figure 9, the natural factors overall have a positive influence on ESs, but the social factors generally have a negative influence on ESs. Dem, Slope and Tem have strong linear characteristics. Slope has the strongest positive impact on ESs, followed by DEM, with Tem exerting the weakest positive influence on ESs. In contrast, Tp and Sand have a strong nonlinear feature. The nonlinear characteristic of Tp gradually intensifies as the independent variable increases. When Tp is below 0.3, the nonlinear effect is not significant. However, as Tp exceeds 0.3, the nonlinear characteristic intensifies, suggesting that high-frequency typhoons have a complex impact on ESs. This effect is particularly strong at higher values of Tp. The nonlinear characteristic of Sand gradually diminishes with the growth of the independent variable. When Sand is below 0.5, its nonlinear features are prominent; however, when Sand exceeds 0.5, these features gradually fade. This phenomenon means that the impact of Sand on ESs is subject to uncertainty when Sand is lower, but sand has a significant negative impact on ESs when sand is higher. Therefore, Tp and Sand exhibit stronger spatial non-stationarity than other natural factors and possess a more complex relationship with ESs.

As far as social factors are concerned, Urban, Pm and Cult exhibit strong linear features and have a negative impact on ESs. Specifically, Cult has a stronger negative influence on ESs, while Urban and Pm have a weaker negative impact on ESs. In contrast, Pop, LST and Port have strong nonlinear characteristics, with Pop and Port showing more pronounced nonlinearity than LST. The nonlinear feature of Pop diminishes as the independent variable increases. When Pop is less than 0.5, its nonlinear characteristic is significant, but as Pop exceeds 0.5, this nonlinearity gradually diminishes. This suggests that lower Pop has a more uncertain impact on ESs, while higher Pop exert a clear inhibitory effect on ESs. The nonlinear feature of LST diminishes as the independent variable increases. LST positively promotes ESs when it is below 0.6 and negatively inhibits ESs when it exceeds 0.6, demonstrating that LST has a dual effect on ESs. Port has the strongest nonlinear property, with its trend showing significant segmentation. When Port ranges from 0.2 to 0.7, it exerts a significant negative influence on ESs; outside this range, it shows varying degrees of positive influence on ESs. Moreover, the relationship between Pop and ESs is relatively complex, and its non-stationarity is stronger than that of other social factors.