Using a Light Gradient-Boosting Machine–Shapley Additive Explanations Model to Evaluate the Correlation Between Urban Blue–Green Space Landscape Spatial Patterns and Carbon Sequestration

Abstract

Global ecosystems are facing challenges posed by warming and excessive carbon emissions. Urban areas significantly contribute to carbon emissions, highlighting the urgent need to improve their ability to sequester carbon. While prior studies have primarily examined the carbon sequestration benefits of single green or blue spaces, the combined impact of urban blue–green spaces (UBGSs) on carbon sequestration remains underexplored. Meanwhile, the rise of machine learning provides new possibilities for assessing this nonlinear relationship. We conducted a study in the Yangzhou urban area, collecting Landsat remote sensing data and net primary productivity (NPP) data at five-year intervals from 2001 to 2021. We applied the LightGBM-SHAP model to systematically analyze the correlation between UBGSs and NPP, extracting key landscape metrics. The results indicated that landscape metrics had varying impacts on NPP. At the patch and type level, the Percentage of Landscape was significantly positively correlated with NPP in green space, while the contiguity index and fractal dimension index favored carbon sequestration under certain conditions. The contribution of blue space was lower, with some indicators exhibiting negative correlations. At the landscape level, the contagion index and aggregation index of UBGS had positive effects on NPP, while the division index and landscape shape index were negatively correlated with NPP. The results enhance the understanding of the relationship between UBGS and carbon sequestration, and provide a reference for urban planning.

1. Introduction

Climate change has emerged as a critical global ecological issue, with carbon emissions being regarded as its main driver [1]. Cities serve as hubs for human economic activities and are significant contributors to fossil energy use and carbon emissions. As the gathering place of human production activities, cities are an important source of fossil energy consumption and carbon emissions [2,3], contributing about 70% of global CO₂ emissions [4], and even up to 80% in China [5]. In its Sixth Assessment Report (AR6), the Intergovernmental Panel on Climate Change (IPCC) emphasized the urgency of cutting urban carbon emissions, in order to achieve the goal of limiting global temperature rise to 1.5 °C by 2025. China, being one of the largest carbon emitters, has pledged to reach a carbon peak by 2030 and attain carbon neutrality by 2060 [6], and has constructed a series of ‘1 + N’ policy systems [7]. Enhancing the carbon sequestration capacity of cities can facilitate offsetting part of the carbon emissions and represents a crucial effort in addressing climate change [8].

Urban Blue–Green Spaces (UBGSs) is a general term for blue spaces and green spaces. Commonly, blue spaces include water bodies such as rivers, lakes, wetlands, and coastal regions; green spaces encompass woodlands, grasslands, croplands, urban parks, and other vegetated areas [9]. With the rapid progress of urbanization, the urbanization rate in China has reached 66.16%, with the greening rate in built-up areas standing at 39.94%, and the area designated for water bodies and water infrastructure covering 36.296 million hectares [10,11,12]. The ecological functions of UBGSs have become increasingly significant [13,14]. Although carbon sequestration functions of large-scale ecosystems (e.g., forests, grasslands, and oceans) have been extensively studied [15,16,17], research on carbon sequestration at the urban scale have primarily focused on individual types of green spaces, often overlooking the synergistic effects of UBGSs [18]. Therefore, it is of great significance to investigate the carbon sequestration benefits of UBGSs.

Gaining a deep understanding of the role of UBGSs in carbon sequestration requires analyzing the factors that influence it. Vegetation types and community structure are widely recognized as key factors affecting carbon storage. For example, Gratani et al. [19] found that the differences in carbon sequestration capacity of urban parks were primarily due to variations in vegetation types; Jo et al. [20] indicated that green spaces with fast growth rates, high density, and diverse structures have stronger carbon storage and absorption capacities. In addition, landscape patterns also play a significant role in carbon sequestration. Ahern [21] pointed out that optimizing the spatial structure and layout of green spaces can enhance carbon sinks and other ecological benefits; Grafius et al. [22] found a positive correlation between patch area of green spaces and carbon sequestration; Godwin et al. [23] revealed a negative correlation between patch edge complexity and carbon density. However, most of the existing studies have focused on single types of green spaces and lack a comprehensive investigation of the synergistic effects of blue–green spaces. In fact, there exists a coupling feedback mechanism between terrestrial carbon absorption and water body availability [24,25], which enhances carbon flux and promotes the carbon cycle [26]. Furthermore, the composite structure of blue–green spaces can further increase carbon sequestration capacity by enhancing biomass accumulation and soil carbon sequestration effects [27,28]. This synergistic effect has been validated in several studies. For example, Jiang et al. [29] found that the importance of waterfront distance factors in carbon sinks reached 28.88%. Li et al. [30] demonstrated that coordinated layouts of blue–green spaces enhance carbon sequestration efficiency. Therefore, a comprehensive study of the synergistic effects of urban blue–green space patterns (UBGSPs) will contribute to a deeper understanding of their role in carbon sequestration.

Carbon sequestration accounting methods mainly include sample plot measurements, model calculations, micrometeorological methods, and remote sensing (RS) estimation methods. Sample plot measurements estimate carbon storage through data collection and formula calculations, but are limited by their inability to provide continuous observation and require destructive experimentation [31]. Model calculations combine plot data with environmental variables (e.g., vegetation characteristics, meteorological data) to assess carbon sequestration. However, common models such as Citygreen, I-Tree Eco, and NTBC require validation when applied to regions with different climatic and physiological conditions [32,33]. Micrometeorological methods estimate carbon flux by measuring CO₂ concentration and meteorological parameters, with Eddy covariance (EC) being a common approach, though it is sensitive to climate and topography [34]. The remote sensing (RS) method estimates regional carbon storage by collecting surface vegetation data via satellites or sensors, offering broader coverage and stronger spatial–temporal continuity compared to ground-based methods, making it ideal for dynamic monitoring of regional carbon sequestration [35,36]. For this study, the RS estimation method was selected to obtain carbon sequestration data at the urban scale.

Although the factors influencing the carbon sequestration capacity of UBGSs have been widely discussed, accurately characterizing the complex relationship between UBGSs and carbon sequestration remains a challenge. The relationship between UBGSs and carbon sequestration is multidimensional and dynamic, influenced not only by temporal and spatial variations [37] but also by environmental threshold effects [38]. Specifically, when certain ecological conditions or land use changes exceed a critical threshold, the carbon sequestration capacity may undergo abrupt shifts or reversals. Therefore, capturing these nonlinear dynamic features is crucial for understanding the carbon sequestration role of UBGS. Traditional methods, such as Pearson correlation coefficient and Spearman’s rank correlation coefficient [39,40], assume linear relationships, often oversimplifying the complex dynamics within ecosystems [41]. In contrast, machine learning methods are able to flexibly capture the relationships between variables, making them better suited for handling complex, large-scale, multidimensional data. In particular, gradient boosting algorithms like LightGBM [42], by capturing complex interactions between variables, can accurately reveal ecological processes, providing more precise modeling methods for understanding carbon sequestration benefits.

While machine learning performs well in complex tasks, it often faces the “black-box” issue [43], which makes it difficult to interpret the model’s predictions. Widely used interpretation methods include the Partial Dependence Plot (PDP), Local Interpretable Model-agnostic Explanations (LIMEs) and Shapley Additive Explanations (SHAP). PDP can demonstrate the impact of individual features or feature combinations, but it is limited in capturing complex feature interactions [44]. LIMEs explain individual predictions through local linear approximations but may overlook global patterns [45]. In contrast, SHAP offers both local and global interpretability by quantifying feature contributions to individual predictions and the overall model [46], making it particularly suitable for understanding complex relationships, such as the one between UBGSs and carbon sequestration.

Yangzhou is a typical water town in China, renowned for its abundant blue–green resources. In recent years, with rapid urban expansion, the carbon emission pressure in Yangzhou has gradually increased [47], posing challenges to the ecosystem. Maintaining ecological balance in Yangzhou and enhancing the city’s carbon sink capacity have become urgent priorities.

This study aims to:

• Explore the complex correlation between UBGSPs and carbon sequestration efficiency.
• Identify key landscape metrics influencing carbon sequestration in UBGSs.
• Propose spatial patterns for high carbon sequestration UBGSs, providing scientific references for optimizing UBGSs configurations and enhancing carbon sequestration potential.

2. Materials and Methods

2.1. Study Area

Located in central Jiangsu Province, Yangzhou (32°15′ N–32°25′ N, 119°01′ E–119°54′ E) spans an area of 6591.21 square kilometers. It is an important city in the Yangtze River Delta urban agglomeration and Nanjing Metropolitan Circle. By the end of 2023, Yangzhou’s urbanization rate had reached 72.79%, with a permanent population of 4.585 million. Yangzhou has a subtropical monsoon climate, marked by mild and humid conditions, with an average annual precipitation of 1020 mm. As a water source for the Eastern Route of the South-to-North Water Diversion Project, the city’s water area accounts for over 25%. The Grand Canal, a UNESCO World Heritage site, runs through the city, earning Yangzhou the title of “the First City of the Grand Canal in China”. The city’s green space coverage reaches 44.9%, with vegetation comprising deciduous broad-leaf forests and evergreen broad-leaf mixed forest belts, forming a favorable blue–green natural landscape.

According to the Territorial Spatial Planning of Yangzhou (2021–2035) [48], the city plans to establish an overall spatial pattern comprising “one district, one belt, one center, and three zones”. “One district” refers to the Yangzhou metropolitan area, “one belt” to the cultural charm belt along the Grand Canal, and “one center” to the ecological green heart of Gaoyou Lake. The “three zones” are composed of distinctive agricultural regions, including the northern Lixiahe region, the southern Yangtze River region, and the western hilly region. We selected the Yangzhou metropolitan area as the research site (Figure 1), characterized by a dense population and frequent human activities. The specific area includes the municipal districts of Yangzhou (Guangling District, Hanjiang District, and Jiangdu District) and Yizheng City, covering a total of 3162.07 square kilometers.

2.2. Data Acquisition and Pretreatment

2.2.1. Remote Sensing

Remote sensing images with a 30 m spatial resolution from the Landsat-5 and Landsat-8 satellites were used [49], covering five periods of multi-band remote sensing data in 2001, 2006, 2011, 2016, and 2021. The combination of Landsat-5 and Landsat-8 ensured data continuity from 2001 to 2021 maintained consistent image quality [50,51]. To balance monitoring requirements and data management needs, a five-year interval was selected to capture long-term trends [52,53]. Next, we used the Seamless Mosaic tool in ENVI 5.6z software to stitch the images from two-track numbers, followed by radiometric calibration and atmospheric geometric correction, to ultimately achieve high-quality remote sensing images.

2.2.2. NPP

The MOD17A3H product from the National Aeronautics and Space Administration (NASA) served as the source for NPP data. This product is based on the Terra/MODIS satellite and provides annual vegetation productivity data at 500 m resolution on a global scale [54]. The dataset is rigorously calibrated for high accuracy and covers multiple years, facilitating long-term time series analyses. It has found extensive application in research fields including ecological monitoring, land use change, and carbon sequestration assessment. We used MRT 4.0 software to process, convert, project, mosaic, and batch process the acquired NPP data. Then, we applied ArcGIS 10.3 to clip the data according to the study area, resulting in the spatial distribution of net primary productivity in the Yangzhou metropolitan area for several specific years between 2001 and 2021 (Figure 2).

2.2.3. Land Use and Blue–Green Spatial Data

The Support Vector Machine (SVM) in the ENVI 5.6 software was used to perform supervised classification on the processed Landsat remote sensing images. SVM is a widely used model for classifying hyperspectral remote sensing images [55], offering high classification accuracy and strong generalization ability, which enables effective distinction of land cover types in high-dimensional data [56]. Moreover, through structural risk minimization, SVM reduces the risk of overfitting, ensuring high reliability and stability in classification tasks [57]. Based on the current land use classification standard (GB/T 21010-2017) [58], land types were classified into five categories: built-up land, forest and grassland, cropland, water bodies, and unused land. The classification results were evaluated using a Confusion Matrix. The Kappa coefficients for the five years ranged from 0.86 to 0.91, meeting the accuracy requirements of the study. Using the reclassification tool in ArcGIS 10.3, water bodies were classified as blue space, while forest, grassland, and cropland land were classified as green space. Built-up land and unused land were categorized as other spaces, resulting in the blue–green space layout of the Yangzhou metropolitan area for the five periods (Figure 3).

2.3. Methods

2.3.1. Landscape Metrics

Landscape metrics are quantitative indicators that condense landscape information [59], and are used to reflect landscape characteristics and their trends. They are applicable to various ecological environment studies, including urban blue–green spaces. Drawing on prior research [22,60,61] and taking into account the unique characteristics of the study area, we selected a total of 12 landscape metrics at three levels: patch, type, and landscape (

Table 1. Blue–green space landscape metrics.

Classification	Metrics	Abbreviations	Formula	Description
Patch level	Dimensionality index	FRAC	$\mathrm{F}\mathrm{R}\mathrm{A}\mathrm{C}={\displaystyle \frac{2\ln \left(0.25{\mathrm{p}}_{\mathrm{i}\mathrm{j}}\right)}{\ln \left({\mathrm{a}}_{\mathrm{i}\mathrm{j}}\right)}}$	Quantifying the shape complexity of patches.
	Related circumscribing circle	CIRCLE	$\mathrm{C}\mathrm{I}\mathrm{R}\mathrm{C}\mathrm{L}\mathrm{E}=1-\left[{\displaystyle \frac{{\mathrm{a}}_{\mathrm{i}\mathrm{j}}}{{\mathrm{a}}_{\mathrm{i}\mathrm{j}}^{\mathrm{s}}}}\right]$	Assessing the degree of circularity of patches.
	Proximity index	CONTIG	$\mathrm{C}\mathrm{O}\mathrm{N}\mathrm{T}\mathrm{I}\mathrm{G}={\displaystyle \frac{\left[\frac{\sum _{\mathrm{r}=1}^{\mathrm{z}}{\mathrm{C}}_{\mathrm{i}\mathrm{j}\mathrm{r}}}{{\mathrm{a}}_{\mathrm{i}\mathrm{j}}}\right]-1}{\mathrm{v}-1}}$	Reflecting the spatial connectivity between patches.
Class level	Area proportion	PLAND	$\mathrm{P}\mathrm{L}\mathrm{A}\mathrm{N}\mathrm{D}={\displaystyle \frac{\sum _{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{i}\mathrm{j}}}{\mathrm{A}}}\times 100$	Patch area percentage, with the largest patch representing the dominant landscape.
	Maximum plaque index	LPI	$\mathrm{L}\mathrm{P}\mathrm{I}={\displaystyle \frac{{\mathrm{m}\mathrm{a}\mathrm{x}}_{\mathrm{j}=1}^{\mathrm{n}}\left({\mathrm{a}}_{\mathrm{i}\mathrm{j}}\right)}{\mathrm{A}}}\times 100$	The proportion of the largest patch to the total landscape area.
	Edge density	ED	$\mathrm{E}\mathrm{D}={\displaystyle \frac{\sum _{\mathrm{k}=1}^{\mathrm{m}}{\mathrm{e}}_{\mathrm{i}\mathrm{k}}}{\mathrm{A}}}\times 10000$	Characterizing the complexity of patch edges.
Landscape level	Separation index	DIVISION	$\mathrm{D}\mathrm{I}\mathrm{V}\mathrm{I}\mathrm{S}\mathrm{I}\mathrm{O}\mathrm{N}=\left[1-{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{m}}}{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{n}}}{\left({\displaystyle \frac{{\mathrm{a}}_{\mathrm{i}\mathrm{j}}}{\mathrm{A}}}\right)}^2\right]$	Indicates the degree of fragmentation of the landscape; a higher value signifies greater fragmentation.
	Perimeter-area fractal dimension	PAFRAC	$\mathrm{P}\mathrm{A}\mathrm{F}\mathrm{R}\mathrm{A}\mathrm{C}={\displaystyle \frac{2/\left\{\mathrm{N}\sum _{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{i}\mathrm{j}}\times \ln {\mathrm{a}}_{\mathrm{i}\mathrm{j}}-\left(\sum _{\mathrm{i}=1}^{\mathrm{m}}\sum _{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{i}\mathrm{j}}\right)\left(\sum _{\mathrm{i}=1}^{\mathrm{m}}\sum _{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{a}}_{\mathrm{i}\mathrm{j}}\right)\right\}}{\left(\mathrm{N}\sum _{\mathrm{i}=1}^{\mathrm{m}}\sum _{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{i}\mathrm{j}}^2\right)-{\left(\sum _{\mathrm{i}=1}^{\mathrm{m}}\sum _{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{i}\mathrm{j}}\right)}^2}}$	Measures the complexity of landscape boundaries and the irregularity of shapes.
	Aggregation index	AI	$\mathrm{A}\mathrm{I}=\left[{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{m}}}\left({\displaystyle \frac{{\mathrm{g}}_{\mathrm{i}\mathrm{i}}}{\mathrm{m}\mathrm{a}\mathrm{x}{\mathrm{g}}_{\mathrm{i}\mathrm{i}}}}\right){\mathrm{p}}_{\mathrm{i}}\right]\times 100$	The degree of landscape patch aggregation, with higher aggregation aiding in maintaining ecological connectivity.
	Landscape Shape Indicators	LSI	$\mathrm{L}\mathrm{S}\mathrm{I}={\displaystyle \frac{\mathrm{E}}{\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{E}}}$	Indicates the complexity of the overall shape of blue–green space; a higher value signifies increased complexity.
	Average patch size	AREA_MN	$\mathrm{A}\mathrm{R}\mathrm{E}\mathrm{A}\_\mathrm{M}\mathrm{N}={\displaystyle \frac{\mathrm{A}}{\mathrm{N}}}{10}^6$	Calculates the average size of patches.
	Contagion index	CONTAG	$\mathrm{C}\mathrm{O}\mathrm{N}\mathrm{T}\mathrm{A}\mathrm{G}=\left\{1+{\displaystyle \frac{\sum _{i=1}^m\sum _{k=1}^m\left[p_i〈\frac{g_{ik}}{\sum _{k=1}^m{g}_{ik}}〉\right]·\left[\mathit{ln}{p}_i〈\frac{g_{ik}}{\sum _{k=1}^m{g}_{ik}}〉\right]}{2\mathit{ln}\left(m\right)}}\right\}\times 100$	Characterizes the aggregation and connectivity of different patches within the landscape.

).

2.3.2. Selection of Optimal Grain Size and the Moving Window Method

The UBGSPs exhibit complexity and multiscale properties, and their impacts on ecological processes or ecosystem services change at different spatial scales [62]. Therefore, selecting the appropriate analysis scale is crucial for accurately capturing how spatial patterns influence ecological processes. To identify an appropriate grain size for UBGSPs, we selected 19 grain sizes for data resampling within the range of 30 m to 150 m and calculated 15 landscape metrics with Fragstats 4.2. By analyzing the curves of the landscape grain size for seven key indices (TA, LPI, SHAPE_MN, COHESION, SPLIT, DIVISION, MESH) [63], we found that most indices stabilized between 45 and 60 m. Combined with the landscape area loss curve (Appendix A), we identified 50 m as the optimal analysis grain size. Based on this, we used the moving window method [64], setting different window radii to calculate landscape metrics. Combined with the semi-variogram [65], we determined the optimal analysis extent. The analysis indicated that a window size of 800 m resulted in minimal spatial variability and strong autocorrelation of the landscape metrics (Appendix A).

2.3.3. LightGBM-SHAP Model

LightGBM (Light Gradient Boosting Machine) is an efficient machine learning algorithm based on the Gradient Boosting Decision Tree (GBDT), developed by Microsoft. It captures complex nonlinear relationships through a progressively weighted decision tree, and shows significant advantages in computational efficiency, accuracy, and noise resistance. LightGBM introduces a histogram-based algorithm [66] that discretizes continuous features into histogram bins for splitting (Figure 4a). This significantly reduces computation time, making it particularly suitable for high-dimensional data and large-scale samples. Unlike traditional level-wise growth strategies, LightGBM adopts a leaf-wise growth strategy (Figure 4b), where it prioritizes expanding leaf nodes that minimize loss the most [67]. This enables the construction of more complex models at the same depth, improving its ability to capture nonlinear features. By combining Gradient-based One-Side Sampling (GOSS) and Exclusive Feature Bundling (EFB) techniques [68], LightGBM enhances computational efficiency and effectively reduces noise, resulting in exceptionally high accuracy and robustness.

This study constructed datasets based on landscape metrics and NPP data from 2001, 2006, 2011, 2016, and 2021. Using ArcGIS 10.3, we created 10,000 random sampling points in green space and 5000 in blue space. We extracted landscape metrics and NPP values at both the patch and class levels for the sampling points, removing outliers to build 10 datasets. In Fragstats, we adjusted the parameters of the classification definition profile to ensure only blue–green spaces were included, creating 15,000 random sampling points and extracting corresponding data at the landscape level to build 5 additional datasets. To guarantee the model’s capacity for generalization, eighty percent of the sampling points were allocated to the training set, while the remaining twenty percent were set aside for the testing set.

The model was optimized through Bayesian opti mization [69]. Key parameters included n_estimators, max_depth, and learning_rate to identify the optimal parameter combination in fewer iterations. Model performance was evaluated using three metrics: R-squared coefficient of determination (R²), mean squared error (MSE), and mean absolute error (MAE), with ten-fold cross-validation used to assess the model’s generalization ability. The formulae for the indicators are as follows:

$${\mathrm{R}}^2=1-{\displaystyle \frac{{\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{{\sum }_{i=1}^n{(y_i-\overline{y})}^2}}$$

(1)

$$\mathrm{M}\mathrm{E}\mathrm{S}={\displaystyle \frac{1}{n}}\sum _{i=1}^n{(y_i-{\widehat{y}}_i)}^2$$

(2)

$$\mathrm{MAE}={\displaystyle \frac{1}{n}}{\sum }_{i=1}^n|y_i-{\widehat{y}}_i|$$

(3)

where $y_i$ is the actual value, ${\widehat{y}}_i$ is the predicted value, $\overline{y}$ is the mean of the actual values, and $n$ is the sample size.

The results showed that both MSE and MAE were low across the 15 datasets, with R² values ranging from 0.856 to 0.922. Results from the ten-fold cross-validation varied between 0.786 and 0.906, demonstrating that the developed LightGBM model met the expected accuracy for both the training and testing sets.

Shapley Additive Explanations (SHAP) is a method for interpreting machine learning model outputs, addressing the “black box” problem. It is based on Shapley values from cooperative game theory [46], assigning importance weights to features and quantifying their contribution to the model’s predictions whether positive or negative. In this study, SHAP was applied to interpret the prediction results of the LightGBM model, aiding in the identification of contributions to carbon sequestration under various landscape features. The formula is as follows:

$${\varphi }_i\left(v\right)={\sum }_{S\subseteq N\backslash \left\{i\right\}}{\displaystyle \frac{\left|S\right|!\left(\left|N\right|-\left|S\right|-1\right)!}{\left|N\right|!}}\left[v\left(S\cup \left\{i\right\}\right)-v\left(S\right)\right]$$

(4)

where ϕi (v) is the Shapley value for feature i; N is the set of all features; Si is the set of variable rankings, and v(s) is the optimal value of the set S.

3. Results

3.1. Spatiotemporal Variation in NPP

From 2001 to 2021, NPP in the study area exhibited a general upward trend (Figure 5a), reflecting an improvement in vegetation growth over the 20-year period. The spatial distribution of NPP was significantly heterogeneous, with high-value areas gradually expanding, particularly in the suburbs and towns in the eastern and northern regions, which were primarily characterized by cropland, forests, and protected areas. The Theil–Sen slope calculation (Figure 5b) and Mann–Kendall test (Figure 5c) [70] revealed that 63.39% of the pixels exhibited a significant increasing trend in NPP, while 3.29% showed a significant decreasing trend. The decreasing pixels were clustered at the boundaries of the central urban area and its sub-districts, indicating that urban expansion negatively influences NPP. Overall, although most regions in the study area experienced an increase in NPP, spatial heterogeneity remained due to urban expansion.

3.2. Variation in Blue–Green Space Spatial Pattern

The findings for UBGSP variation in the study area (Appendix B) indicate that, at the patch level, changes in indices for green space and blue space were not significant at the patch level. At the class level, the PLAND, LPI, and ED of green space declined overall, reflecting reduced fragmentation and alleviated edge effects. In contrast, three indices for blue space slightly increased, indicating improvements in area and connectivity, along with enhanced edge effects. At the landscape level, AREA_MN increased, suggesting an overall reduction in fragmentation. The increasing and then decreasing trend of LSI and PAFRAC indicated a decrease in the degree of discreteness. CONTAG showed a decreasing trend, suggesting a weakening of patch connectivity in the UBGSs. Overall, through the integration of urban green spaces and water bodies, Yangzhou reduced landscape fragmentation and improved the connectivity and aggregation of blue–green spaces.

3.3. The Correlation Between UBGSs and Carbon Sequestration Based on LightGBM-SHAP

3.3.1. Global Interpretation

The results of the LightGBM-SHAP model were as follows. The SHAP summary plot (Figure 6) was used for global interpretation of the model, displaying the importance ranking of features within the dataset. The ranking indicated importance from high to low.

In green space, the top three important indices from 2006 to 2016 were PLAND, CONTIG, and FRAC, while in 2001 and 2021, they were PLAND, CONTIG, and ED. In general, PLAND had the greatest influence on NPP, demonstrating a clear positive correlation, highlighting the beneficial impact of large green patches on NPP. CONTIG ranked second in contribution, also exerting a positive effect on NPP. The SHAP value rankings for FRAC and ED indicated that both of them had a strong influence on NPP across different years.

In the blue space, the top three contributing indices in 2001, 2011, and 2016 were CONTIG, PLAND, and LPI, while in 2006, they were PLAND, CONTIG, and ED, and in 2021, CONTIG, PLAND, and FRAC. Overall, the CONTIG index exerted the strongest influence on NPP, with PLAND and LPI following in impact.

The landscape-level indices of blue–green space displayed a complex relationship with NPP. In the SHAP summary plot from 2001 to 2021, the three most important contributing indices were DIVISION, LSI, and CONTAG, with LSI declining in 2021, when AREA_MN became the third most important index. This indicated that the fragmentation, shape complexity, and connectivity of UBGS were key factors influencing NPP.

3.3.2. Local Interpretation

The SHAP dependence plot provided a local interpretation, illustrating how individual features influenced the prediction results and revealing interactions between features. We selected the metrics with higher contribution rankings from each year in the summary plot for analysis.

(1)

Patch and class level analyses

In green space (Figure 7), a significant positive relationship was observed between PLAND and NPP, indicating a beneficial effect. The trend of CONTIG was not significant when its value was low, but a clear positive trend was observed within the range of [0.7, 1]. The FRAC values took 1.2 as the critical point, with opposite trends on either side affecting NPP. When FRAC < 1.2, the SHAP value gradually decreased as FRAC increased; conversely, when FRAC > 1.2, the SHAP value increased, moving towards positive values. The trend of ED was more significant in 2001 and 2021; as ED increased, the SHAP value decreased, indicating a negative impact on NPP.

Compared to green space, the overall effect of blue space on NPP was weaker (Figure 8). CONTIG contributed positively to NPP in the low to medium value range; however, when it exceeded 0.9, it had a significant negative impact. Although some PLAND values were dispersed, there was an overall negative linear trend with NPP. Low LPI values clustered around a SHAP value of 0, and as the LPI increased, its SHAP value exhibited an upward trend, though this was not significant. When the ED value was below 40, most SHAP values were greater than 0. With an increase in the ED value, the SHAP value exhibited a downward trend.

At the patch and class level, the relationship between landscape metrics and NPP was stronger in green space than in blue space. Given the complex interactions between green and blue spaces, their effects on NPP are not isolated. Therefore, taking a landscape-level perspective is crucial for a more comprehensive understanding of how combinations of UBGSs impact carbon sequestration.

(2)

Landscape level analyses

The SHAP dependence plot analysis revealed the contribution mechanisms of UBGSP to NPP (Figure 9), which we analyzed according to the comprehensive ranking of contribution from highest to lowest in the summary plot.

The DIVISION value represented the separation degree of patches, ranging from 0 to 1. This range indicated a transition from a single green or blue patch to an increasing number of both types, eventually becoming fragmented by other patch types. DIVISION showed a significant negative correlation with SHAP values, suggesting a negative effect on NPP. LSI measured the complexity of patch shapes and was typically associated with landscape fragmentation. It influenced NPP similarly to DIVISION; when the value was around three, it served as a boundary point for SHAP values. Higher LSI values corresponded to lower SHAP values, demonstrating a negative correlation. The trend of CONTAG fluctuated over the years, with a more marked positive effect on NPP during the notable years of 2001, 2011, and 2021. AI was used to measure the aggregation degree of a particular type of patch (green space or blue space), and it showed a typical threshold effect. When the AI value was in the range of [71,72], the SHAP value was relatively high and positive. As the AI values deviated from this range, SHAP values decreased. AREA_AM was positively correlated with NPP overall, with higher AREA_AM values typically corresponding to higher NPP values. In addition, the contribution of PAFRAC to NPP was low, with no clear trend in its effect.

In summary, the connectivity and aggregation of blue–green spaces positively influenced NPP, while fragmentation and high separation brought negative impacts. Compared to isolated green or blue spaces, the synergistic interaction of UBGS produced more favorable effects. In the integrated ecological environment where green and blue spaces coexist, a more conducive habitat for plant growth is created, thereby enhancing carbon sequestration capacity. This highlights the importance of optimizing the layout and comprehensive management of UBGSs.

4. Discussion

4.1. Impact of Green Space and Blue Space Spatial Patterns on Carbon Sequestration

Numerous studies have confirmed the positive influence of green space on NPP. We found that PLAND, CONTIG, FRAC, and LPI in green spaces exhibited positive correlations with NPP, while PLAND and CONTIG having the strongest impacts. This indicates that larger patch sizes and higher connectivity can greatly improve carbon sequestration capacity [61,73]. Large green patches provide plants with more sunlight, water, and nutrients, creating stable habitats, while high connectivity promotes material flow within the ecosystem. Together, these factors greatly enhance NPP. On the other hand, the results for FRAC indicate that higher patch complexity contributes to enhanced carbon sequestration capacity. This is primarily because complex patch boundaries are more irregular, creating diverse habitats that facilitate vegetation growth and carbon storage. In contrast, ED was negatively correlated with NPP, indicating that edge effects may hinder carbon sequestration. This supports the results of Meeussen et al. [74], who observed that edge effects can have negative impacts on carbon sequestration under specific conditions. Due to the swift urban growth in recent years, green spaces located on the outskirts of Yangzhou’s urban areas, despite having long interfaces, have been more susceptible to human degradation, leading to suboptimal NPP outcomes.

Our findings on the relationship between blue space connectivity and NPP differ from previous studies [75]. The analysis revealed that highly connected water systems do not always correlate with higher NPP. Although the Grand Canal in Yangzhou is well-connected, human disturbances have affected it. The hardened banks hinder the exchange between blue and green spaces, limiting vegetation growth and resulting in habitat simplification, which lowers NPP values. The findings align with those of Wohl et al. [76], who noted that human activities disrupt natural flow patterns, weakening the organic carbon storage capacity of river corridors. This also explains why edge benefits are negatively correlated with NPP. Additionally, the negative correlation between PLAND and NPP may be attributed to the dominant role of green plants in carbon sequestration, while the amount of vegetation in blue space is limited. Even with an expansion of water areas, the increase in aquatic plant biomass is minimal, which could lead to insufficient NPP growth.

4.2. Impact of Blue–Green Space Spatial Patterns on Carbon Sequestration

The impact of spatial patterns in combined UBGSs on NPP has been less explored in the existing literature, our study offers new insights into this area. The results for UBGSs confirmed the positive impact of their synergy on carbon sequestration. DIVISION was negatively correlated with NPP, indicating that fragmentation of UBGSs reduces carbon sequestration capacity. Higher DIVISION typically corresponds to smaller patches, which are more vulnerable to the island effect, such as reduced species diversity and weakened ecological functions. The northern area of Yangzhou Biodiverse and Sci-Tech City once had large patches of wetlands and farmland. However, with urbanization and the development of infrastructure like roads, blue–green spaces were fragmented into smaller patches, resulting in a decrease in NPP in this region. On the other hand, when the AI value was in the middle range, the aggregation of blue–green spaces was highest, and NPP peaked. In these spaces, the intermixing of water bodies, wetlands, and vegetation expanded the core areas, minimizing edge effects. This promotes a stable ecological environment and increases biodiversity, ultimately enhancing carbon sequestration by vegetation [71,77]. Slender West Lake in Yangzhou serves as a typical example. In addition, the negative correlation of LSI with NPP is in agreement with the results reported by Yuan et al. [60] in Nanjing. Complex, irregularly shaped patches have longer edges, making them more susceptible to external factors, which weaken the stability of internal habitats. This explains the lower NPP observed in Yangzhou’s urban expansion areas. Furthermore, the high connectivity of blue–green spaces contributes to enhanced carbon sequestration, aligning with previous research in ecosystems such as forests, grasslands, and rivers [78,79]. For example, along the Yangtze River, certain areas have maintained a natural ecological environment with good connectivity. This promotes the flow of ecological materials, improving resource utilization efficiency and thus supporting the carbon sequestration effect of vegetation. AREA_AM and NPP showed a positive connection, suggesting that larger blue–green space patches were more beneficial for NPP accumulation. We believe that larger patches provide more extensive habitats and richer resources, while also reducing the impact of edge effects and enhancing ecological stability, all of which contribute to vegetation’s capacity for carbon sequestration [80]. The limited contribution of PAFRAC may be due to the fact that the carbon sequestration capacity of the UBGSs relies more on the combined effects of terrestrial and aquatic ecosystems, rather than the shape of the patches being a key factor.

4.3. Optimization Strategies for UBGSP

Achieving “carbon neutrality” and a “carbon peak” has become a strategic goal for nations in addressing global warming, with UBGSs playing a key role in this effort. Using Yangzhou’s 2021 blue–green spatial layout as an example, we extracted typical high-carbon sequestration UBGS units with a 5 km × 5 km window (

Table 2. Typical high-carbon sequestration UBGS units in study area.

UBGS Characteristics	5 km × 5 km Unit	Green Space Patch and Class-Level Metrics	Blue Space Patch and Class-Level Metrics	Blue–Green Space LANDSCAPE-Level Metrics
Large-scale UBGSs	West Suburb Forest Park	1.00 < FRAC < 1.16 0.36 < CIRCLE < 0.68 0.41 < CONTIG < 0.8 PLAND = 75.47 LPI = 75.35 ED = 70.67	1.00 < FRAC < 1.11 0.36 < CIRCLE < 0.71 0.41 < CONTIG < 0.49 PLAND = 1.28 LPI = 0.16 ED = 6.27	LSI = 10.42 AREA_MN = 8.40 PAFRAC = 1.47 CONTAG = 50.69 DIVISION = 0.43 AI = 81.40
High-density UBGSs	Shaobo Lake	1.00 < FRAC < 1.29 0.43 < CIRCLE < 0.79 0.17 < CONTIG < 0.85 PLAND = 50.19 LPI = 33.64 ED = 68.25	1.00 < FRAC < 1.22 0.49 < CIRCLE < 0.83 0.13 < CONTIG < 0.91 PLAND = 24.12 LPI = 10.65 ED = 34.85	LSI = 13.12 AREA_MN = 9.30 PAFRAC = 1.51 CONTAG = 23.55 DIVISION = 0.86 AI = 76.17
Highly aggregated UBGSs	Slender West Lake	1.00 < FRAC < 1.24 0.36 < CIRCLE < 0.85 0.08 < CONTIG < 0.82 PLAND = 25.05 LPI = 8.50 ED = 53.77	1.00 < FRAC < 1.17 0.56 < CIRCLE < 0.90 0.27 < CONTIG < 0.74 PLAND = 3.92 LPI = 2.21 ED = 8.65	LSI = 8.70 AREA_MN = 2.78 PAFRAC = 1.44 CONTAG = 50.20 DIVISION = 0.47 AI = 84.93

) and analyzed their spatial pattern characteristics. Based on these findings, optimization strategies for UBGSs are proposed.

(1)

Enhancing the aggregation and connectivity of UBGSs

The degree of aggregation and connectivity of UBGSs is crucial for ecosystem stability and carbon sink capacity. During urban development, emphasis should be placed on maintaining extensive, concentrated blue and green spaces. The existing key areas, such as parks and wetlands, should be protected from the interference of non-ecological land expansion by designating ecological protection areas [81]. In addition, the urban construction process should focus on the construction of ecological corridors [72]. By integrating urban greenways and rivers, ecological patches such as scattered green spaces, wetlands, and water bodies can be organically connected. Ecological tunnels should be set up in areas separated by roads and bridges to avoid the division of UBGSs by infrastructure. This will result in blue–green spaces with stronger carbon sequestration capacity.

(2)

Optimizing the shape and structure of UBGSs

To improve the carbon sequestration capacity of UBGSs, natural shorelines should be protected and restored in areas along lakes and rivers, preserving natural slopes with vegetation cover [67]. In areas requiring flood control, the use of hard engineering can be reduced through resilient landscape design. The design of UBGSs should focus on a staggered land–water layout, creating moderate zigzag transition zones, to enhance ecological functions. Care should be taken to avoid overly complex shoreline designs, as it may lead to space fragmentation. Additionally, UBGS boundaries should be delineated as ecological buffer zones with multi-layered vegetation, such as grass, shrubs, and trees. This action helps shield the core area from the harmful effects of edge effects and human activities.

(3)

Establishing a multi-scale integrated UBGS network

At the macroscopic level, a continuous urban ecological network [82] can be created by integrating large-scale ecological corridors, green spaces, and water systems, to effectively connect the city with surrounding natural areas. At the meso level, the layout of UBGSs between urban functional areas should be optimized to reduce ecological interference between different land use types. At the micro level, local facilities such as community parks, green roofs [83], and vertical greening can connect scattered ecological patches, enhancing the clustering and ecological value of blue–green spaces. This multi-scale strategy not only increases the carbon Sequestration potential of cities, but also enhances their ecological resilience and provide more social service functions.

4.4. Limitations and Future Research Directions

This study has some limitations. First, the relatively low spatial resolution of MODIS limits the accurate capture of small-scale landscape changes, which could affect the precision of carbon sequestration estimates. Second, while the analytical framework of this study primarily focused on the relationship between UBGSPs and carbon sequestration, carbon sequestration is influenced by other factors such as climate conditions, soil properties, plant species, and vegetative, which were not fully considered in this framework. Lastly, although this study provides theoretical support for research on carbon sequestration in UBGSs, translating the findings into specific policy recommendations and planning guidelines remains challenging.

As technology continues to advance, future research could improve the accuracy and reliability of carbon sequestration estimates by utilizing higher-resolution remote sensing data, combining ecosystem carbon simulation models, and employing multi-source data [84]. Additionally, by incorporating more environmental variables, a comprehensive carbon sequestration evaluation framework could be developed, enabling integrated analysis across spatial and temporal scales. Moreover, future studies should also focus on the translation of research into practical applications, establishing coordination mechanisms among urban management departments (e.g., urban construction, water resources, agriculture, forestry, etc.) to facilitate the integration of carbon sequestration insights into urban policy and planning processes.

5. Conclusions

This study focused on the metropolitan area of Yangzhou and employed the LightGBM-SHAP model to explore the relationship between UBGS and carbon sequestration. Based on typical high-carbon sequestration units, optimization strategies for the UBGS layout in Yangzhou were proposed. The main findings are as follows:

(1)

The relationship between UBGSPs and carbon sequestration was nonlinear, with different landscape metrics demonstrating varying effects at different levels. At the patch and class levels, PLAND, CONTIG, FRAC, and LPI had a more significant effect on carbon sequestration; at the landscape level, DIVISION, LSI, and CONTAG played an important role in carbon sequestration.

(2)

Compared to single green or blue spaces, the synergistic effects of blue–green spaces had a more positive impact on carbon sequestration. UBGSPs with high connectivity (CONTAG) and large areas (AREA_MN) were positively correlated with NPP, while the aggregation index (AI) specifically exhibited a significant threshold effect, with NPP reaching its maximum within the [71,72] range. Conversely, higher values of DIVISION and LSI (>3) were negatively correlated with NPP, indicating that fragmented and irregular blue–green space patches were unfavorable for carbon sequestration.

(3)

Large-scale, high-density, and highly aggregated UBGSs were typical high-carbon sink areas. Under the premise of ensuring appropriate aggregation and connectivity, optimizing patch shape and improving the connectivity and stability of the ecological network helped enhance carbon sequestration capacity of UBGSs.

This study fills the gap in urban carbon sequestration research by incorporating the synergies of UBGS, highlighting their positive impact on carbon sequestration. This finding contributes to carbon sequestration theory, advances urban ecosystem research, and provides a new framework for exploring the interactions between UBGS. Furthermore, the results offer valuable scientific support for urban planning and ecological development, promoting green, low-carbon urbanization, especially in rapidly urbanizing countries like China.