Identification of Urban Green Space Types and Estimation of Above-Ground Biomass Using Sentinel-1 and Sentinel-2 Data

Abstract

High-quality urban green space supports the healthy functioning of urban ecosystems. This study aimed to rapidly assess the distribution, and accurately estimate the above-ground biomass, of urban green space using remote sensing methods, thus providing a better understanding of the urban ecological environment in Xuzhou for more effective management. We performed urban green space classifications and compared the performance of Sentinel-2 MSI data and Sentinel-1 SAR data and combinations, for estimating above-ground biomass, using field data from Xuzhou, China. The results showed the following: (1) incorporating an object-oriented method and random forest algorithm to extract urban green space information was effective; (2) compared with stepwise regression models with single-source data, biomass estimation models based on multi-source data provide higher estimation accuracy (R² = 0.77 for coniferous forest, R² = 0.76 for shrub-grass vegetation, R² = 0.75 for broadleaf forest); and (3) from 2016 to 2021, urban green space coverage in Xuzhou decreased, while the total above-ground biomass increased, with higher average above-ground biomass in broadleaf forests (133.71 tons/ha) compared to coniferous forests (92.13 tons/ha) and shrub-grass vegetation (21.65 tons/ha). Our study provides an example of automated classification and above-ground biomass mapping for urban green space using multi-source data and facilitates urban eco-management.

1. Introduction

In the urgent context of global warming, cities, with their intensive energy consumption and greenhouse gas emissions [1,2], are required to take responsibility for carbon dioxide emission reduction [3,4,5,6]. Urban green spaces (UGSs) are woody and herbaceous vegetated areas open to the public and have valuable natural carbon sequestration functions, as well as playing an important role in the local carbon balance in urban areas [7,8,9,10]. One of the important indicators for quantifying the carbon sequestration capacity of UGSs is above-ground biomass (AGB) [11,12,13]. Rapid mapping of UGSs and the refined acquisition of AGB data for continuous monitoring of green space vegetation are needed, as these will improve scientific exploration of urban ecosystem material cycles and green space carbon sequestration potential. They also facilitate in-depth consideration about UGS management.

Traditionally, accurate information on the distribution of UGS was obtained through field surveys at a considerable expense in terms of labor and time [14,15]. Automatic information acquisition and repeatable observations can be used to monitor UGS, by employing remote sensing image classification technology based on pixels or objects [16,17,18]. In addition, the object-oriented method focuses on the relationship between pixels and their neighbors, which can provide higher accuracy in UGS classification results, as an alternative to the traditional pixel-based image classification method, that can result in the “salt and pepper effect” [19,20,21]. Likewise, it is inappropriate, as well as ecologically destructive, to rely on the traditional vegetation harvesting method, i.e., by harvesting-drying-weighing, for AGB monitoring in environmentally complex and easily disturbed UGSs [22]. The remote sensing method facilitates AGB monitoring without destroying any existing vegetation, while striking a good balance between temporal and spatial scales [23,24]. Its basic idea is to establish a correlation between remote sensing image features and sample point-based measured data, and then use parametric models or machine learning methods to calculate AGB [25,26,27]. Multiple linear regression is an important method for inferring vegetation AGB because it is highly interpretable and understandable [28,29].

Compared to the widely used Landsat satellite data, Sentinel-2 MSI data, freely available from the European Space Agency (ESA) since June 2015, has a higher spatial resolution of 10 m. This improves the detection of vegetation-soil mixed pixels and makes UGS highly separable from other urban surface objects, so it has increasingly become one of the most valuable data sources [30]. In addition, Sentinel-2 MSI data uniquely contains three bands within the red-edge range, which improves retrieval of vegetation canopy information [31,32]. However, the Sentinel-1 radar system, equipped with a C-band (wavelength ~6 cm) sensor, is superior for detecting ground objects without the hindrance of clouds and fog [33]. It reflects vertical structure information and vegetation function in urban green ecosystems [34], which somewhat compensates for the lack of biomass information about branches and stems in the vertical direction of vegetation in the Sentinel-2 MSI data [35]. Even though optical remote sensing data and radar data are the two main data sources for vegetation remote sensing applications, so far remarkable success has been achieved in monitoring natural forests [36,37,38] and grasslands [39,40], and limited emphasis has been placed on smaller but complicated UGS areas [41]. At the same time, it has been used much less frequently than optical remote sensing images, most likely because of the complexity of radar data processing procedures [28]. As such, the potential for combining high-resolution Sentinel-2 and Sentinel-1 multi-source remote sensing data for green space monitoring in urban areas remains to be explored.

The purpose of this study was to explore the potential of integrating Sentinel-2 MSI data and Sentinel-1 SAR data in UGS classification and AGB estimation. We first created Xuzhou UGS distribution maps for three time periods (2016, 2018, and 2021) using the object-oriented classification method. Subsequently, we sought to quantify AGB in UGS by vegetation types. The results of this study improve the understanding of UGS distribution and dynamic changes of AGB in Xuzhou, which will contribute to analyzing the urban carbon budget and provide a scientific basis for future UGS construction.

2. Materials and Methods

2.1. Study Area and Field Data Collection

Known as the capital of China’s construction machinery and the second-largest railway hub in China, Xuzhou (116°22′~118°40′ E, 33°43′~34°58′ N) (Figure 1e) is an important gateway city to the northwest Jiangsu Province, a province bordering the Yellow Sea in eastern China. Influenced by the southeast monsoon, it features a warm temperate semi-humid monsoon climate with sufficient sunlight and moderate rainfall. The annual average temperature and precipitation in Xuzhou are 14.50 °C and 841.20 mm, respectively. As an old resource-based city, Xuzhou has been committed to the construction of UGS in recent years. It is worth mentioning that its urban forest coverage rate ranks among the best in Jiangsu Province, reaching 27.80% in 2021 [42].

The study area (117°2′24.58″~117°19′57.11″ E, 34°8′34.38″~34°23′4.10″ N) is located within the outer ring expressway around Xuzhou (represented by the red line in Figure 1d), which is the living and working circle of urban residents in the central downtown. It covers a geographical area of approximately 527.63 km². The dominant forest species in the urban area is warm temperate deciduous broadleaf [43], and popular tree species are Cinnamomum camphora (L.) Presl., Ginkgo biloba L., and Koelreuteria paniculata Laxm. The endemic herbaceous plants include Ferula licentiana Hand.-Mazz. var. tunshanica (Su) Shan et Q.X.Liu, Rehmannia glutinosa (Gaert.) Libosch.ex Fisch.et Mey.

From October to early December 2017, we conducted field surveys on typical green spaces in Xuzhou. To match the spatial resolution of the Sentinel-2 MSI data, rectangular green quadrats were set up with a fixed size of 10 m × 10 m. The growth indicators of the vegetation in the quadrats were measured and recorded, including diameter at breast height (DBH), height (from the base to the crown), and crown width of each tree; base diameter, height, and crown width (both south-north and east-west directions) of each shrub; and grass coverage and plant height. Afterward, a dataset of AGB with a total of 140 quadrats was created. The hierarchical summary method was used to calculate AGB for trees, shrubs, and grass layers separately, and then they were summed to obtain the AGB in the same quadrat. Furthermore, the AGB of trees and shrubs was estimated indirectly using allometric models (see

Table A1. Tree allometric growth equation for calculating quadrat above-ground biomass of trees.

Tree Species	Biomass Model	R2	References
Platycladus orientalis (L.) Franco	${\mathrm{W}}_{\mathrm{S}}=0.5730{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.8657}$	0.97	[75]
	${\mathrm{W}}_{\mathrm{B}}=0.0043{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.1085}$	0.89
	${\mathrm{W}}_{\mathrm{L}}=0.0038{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.0385}$	0.84
Metasequoia glyptostroboides Hu & W.C. Cheng	${\mathrm{W}}_{\mathrm{S}}=-0.656+0.028{\mathrm{D}}^2\mathrm{H}$	/	[76]
	${\mathrm{W}}_{\mathrm{B}}=-1.258+0.007{\mathrm{D}}^2\mathrm{H}$	/
	${\mathrm{W}}_{\mathrm{L}}=-0.004+0.001{\mathrm{D}}^2\mathrm{H}$	/
	${\mathrm{W}}_{\mathrm{P}}=-0.135+0.003{\mathrm{D}}^2\mathrm{H}$	/
Cupressus funebris Endl.	${\mathrm{W}}_{\mathrm{T}}=0.02479{\mathrm{D}}^{2.0333}$	/	[76]
Pinus	$\ln {\mathrm{W}}_{\mathrm{S}}=-1.710+0.688\ln \left({\mathrm{D}}^2\mathrm{H}\right)$	/	[77]
	$\ln {\mathrm{W}}_{\mathrm{B}}=0.089+0.096\ln \left({\mathrm{D}}^2\mathrm{H}\right)$	/
	$\ln {\mathrm{W}}_{\mathrm{L}}=0.335+0.006\ln \left({\mathrm{D}}^2\mathrm{H}\right)$	/
Magnolia grandiflora L.	${\mathrm{W}}_{\mathrm{S}}=0.0649{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.8131}$	0.97	[78]
	${\mathrm{W}}_{\mathrm{B}}=0.0431{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.6697}$	0.90
	${\mathrm{W}}_{\mathrm{L}}=0.0254{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.8701}$	0.84
Populus L.	${\mathrm{W}}_{\mathrm{S}}=0.0074046{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.0690}$	/	[76]
	${\mathrm{W}}_{\mathrm{B}}=0.0041773{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.9911}$	/
	${\mathrm{W}}_{\mathrm{L}}=0.0715320{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.4489}$	/
Cinnamomum camphora (L.) Presl.	${\mathrm{W}}_{\mathrm{S}}=0.0914{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.7755}$	0.99	[78]
	${\mathrm{W}}_{\mathrm{B}}=0.0099{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.0256}$	0.95
	${\mathrm{W}}_{\mathrm{L}}=0.0011{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.1713}$	0.94
Ginkgo biloba L.	$\ln {\mathrm{W}}_{\mathrm{S}}=-2.23+2.17\ln \mathrm{D}$	0.99	[79]
	$\ln {\mathrm{W}}_{\mathrm{B}}=-6.95+3.36\ln \mathrm{D}$	0.87
	$\ln {\mathrm{W}}_{\mathrm{L}}=-5.21+2.37\ln \mathrm{D}$	0.87
Koelreuteria paniculata Laxm.	${\mathrm{W}}_{\mathrm{S}}=0.587+0.006{\mathrm{D}}^2\mathrm{H}$	0.63	[80]
	${\mathrm{W}}_{\mathrm{B}}=0.059+0.001{\mathrm{D}}^2\mathrm{H}$	0.92
	${\mathrm{W}}_{\mathrm{L}}=0.029+0.001{\mathrm{D}}^2\mathrm{H}$	0.95
$\mathrm{Ligustrum}\mathrm{lucid}$um Ait.	${\mathrm{W}}_{\mathrm{S}}=0.443{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.608}$	0.99	[81]
	${\mathrm{W}}_{\mathrm{B}}=0.287{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.548}$	0.96
	${\mathrm{W}}_{\mathrm{L}}=0.844{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.268}$	0.97
Other hard broadleaf trees	${\mathrm{W}}_{\mathrm{T}}=0.07112{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.910358078}$	/	[76]
Other soft broadleaf trees	${\mathrm{W}}_{\mathrm{S}}=0.012540{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.144}$	/	[76]
	${\mathrm{W}}_{\mathrm{B}}=0.004786{\left({\mathrm{D}}^2\mathrm{H}\right)}^{1.006}$	/
	${\mathrm{W}}_{\mathrm{L}}=0.047180{\left({\mathrm{D}}^2\mathrm{H}\right)}^{0.769}$	/

Note: D represents diameter at breast height; H represents height; ${\mathrm{W}}_{\mathrm{S}}$, ${\mathrm{W}}_{\mathrm{B}}$, and ${\mathrm{W}}_{\mathrm{L}}$ refer to organ component biomass of stems, branches, and leaves; and ${\mathrm{W}}_{\mathrm{T}}$ refers to total above-ground biomass.

and

Table A2. Shrub allometric growth equation for calculating quadrat above-ground biomass of shrubs.

Shrub Species	Biomass Model	R2	References
Ligustrum quihoui Carr.	${\mathrm{W}}_{\mathrm{B}}=26.332 {(\mathrm{C}\mathrm{H})}^{0.666}$	0.95	[82]
	${\mathrm{W}}_{\mathrm{L}}=14.640{\mathrm{C}}^{1.164}$	0.97
Nandina domestica Thunb.	${\mathrm{W}}_{\mathrm{B}}=75.700 {(\mathrm{C}\mathrm{H})}^{1.110}$	0.98	[83]
	${\mathrm{W}}_{\mathrm{L}}=11.109+17.911\ln \mathrm{H}$	0.97
Forsythia suspensa (Thunb.) Vahl	${\mathrm{W}}_{\mathrm{B}}=0.385 {({\mathrm{D}}^2\mathrm{H})}^{1.025}$	0.99	[83]
	${\mathrm{W}}_{\mathrm{L}}=0.187 {({\mathrm{D}}^2\mathrm{H})}^{0.868}$	0.99
Lagerstroemia indica L.	${\mathrm{W}}_{\mathrm{B}}=30.213{\mathrm{H}}^{6.318}$	0.99	[83]
	${\mathrm{W}}_{\mathrm{L}}=6.656{\mathrm{H}}^{5.065}$	0.99
Photinia serratifolia (Desfontaines) Kalkman	${\mathrm{W}}_{\mathrm{B}}=0.310 {({\mathrm{D}}^2\mathrm{H})}^{1.097}$	0.99	[83]
	${\mathrm{W}}_{\mathrm{L}}=0.264 {({\mathrm{D}}^2\mathrm{H})}^{0.916}$	0.99
Tree-like shrub	${\mathrm{W}}_{\mathrm{T}}=0.182{\mathrm{D}}^{2.487}$	0.97	[84]
Other shrubs	${\mathrm{W}}_{\mathrm{T}}=0.429{\mathrm{D}}^{1.386}{\mathrm{H}}^{0.978}$	0.86	[85]

Note: D represents basal diameter; H represents height; C represents crown width (being a mean of both south-north and east-west crown diameter; ${\mathrm{W}}_{\mathrm{B}}$, ${\mathrm{W}}_{\mathrm{L}}$, refers to organ component biomass of branches and leaves; ${\mathrm{W}}_{\mathrm{T}}$ Refers to total above-ground biomass.

, Appendix A), which were established based on correlations between the dry weight of each organ component and easily measurable variables (e.g., DBH, height) for trees [44,45]. The AGB for bamboo forest was also calculated to be 22.50 kg per plant [46], and 8 tons/ha was taken as a reference for calculating the measured grassland biomass in Xuzhou [47]. The biomass calculation results from field data are listed in

Table A3. The above-ground biomass measurements for 140 quadrants.

Id	AGB (kg)	Id	AGB (kg)	Id	AGB (kg)	Id	AGB (kg)	Id	AGB (kg)
1	8.00	2	8.00	3	8.00	4	8.00	5	8.00
6	8.00	7	8.00	8	8.00	9	8.00	10	8.00
11	91.60	12	94.08	13	96.75	14	117.76	15	127.91
16	133.23	17	160.61	18	175.25	19	186.37	20	202.69
21	214.39	22	230.27	23	237.78	24	254.89	25	281.64
26	336.18	27	352.39	28	357.28	29	368.60	30	375.56
31	386.50	32	417.02	33	473.38	34	506.82	35	520.87
36	531.62	37	533.54	38	558.13	39	561.18	40	564.04
41	585.47	42	613.15	43	632.29	44	641.06	45	656.58
46	673.60	47	679.06	48	679.67	49	682.36	50	722.11
51	723.64	52	729.26	53	729.85	54	734.95	55	745.40
56	759.98	57	772.86	58	808.19	59	820.15	60	821.69
61	853.74	62	885.80	63	902.78	64	915.33	65	915.70
66	930.44	67	940.77	68	949.47	69	950.03	70	950.84
71	950.94	72	953.26	73	953.39	74	969.34	75	970.26
76	983.34	77	989.02	78	993.63	79	997.36	80	1001.83
81	1007.31	82	1034.85	83	1037.91	84	1049.60	85	1066.50
86	1071.04	87	1087.98	88	1092.26	89	1106.63	90	1115.03
91	1152.97	92	1153.15	93	1166.04	94	1208.30	95	1210.57
96	1225.00	97	1240.27	98	1247.29	99	1267.29	100	1276.68
101	1285.56	102	1314.06	103	1316.19	104	1332.54	105	1354.43
106	1355.76	107	1380.25	108	1390.47	109	1401.85	110	1412.17
111	1425.54	112	1468.00	113	1477.88	114	1559.17	115	1584.37
116	1594.27	117	1623.34	118	1711.02	119	1735.48	120	1748.48
121	1760.95	122	1806.12	123	1811.93	124	1843.84	125	1848.14
126	1948.03	127	1948.30	128	2013.30	129	2106.08	130	2208.01
131	2217.66	132	2346.12	133	2360.55	134	2379.07	135	2588.19
136	2923.15	137	3296.45	138	3399.43	139	10,331.06	140	11,250.00

, Appendix A.

2.2. Satellite Data

The Sentinel-2 MSI data used in this study were all acquired during the summer months (from June to August), except for images used for modeling that coincided with the autumn field survey period (

Table 1. Remote sensing images used in this study.

Mission	Image Acquisition Date	Image Product Grade	Pixel Spacing (m)
Sentinel-2A/B	28 August 2016	MSI_Level-1C	Bands 2, 3, 4, 8: 10. Bands 5, 6, 7, 8a, 11, 12: 20. Bands 1, 9, 10: 60.
	27 October 2017
	14 June 2018
	23 June 2021
Sentinel-1A	29 August 2016	IW_GRDH	10
	23 October 2017
	8 June 2018
	28 June 2021
SRTM	23 September 2014	SRTM1	30

). Sentinel-2 MSI data pre-processing was mainly performed using Sen2cor 02.05.05 for atmospheric correction and then the data were resampled to 10 m in the SNAP 8.0 (Sentinel Application Platform) environment. Excluding three bands (1, 9, and 10) that had only marginal relevance to vegetation monitoring, the remaining 10 single bands were combined in ENVI 5.3 followed by masking and image mosaicking. We also used the Sentinel-1 SAR data with the dual-polarization VV + VH mode. This product had already been multi-looked and projected to ground range using an Earth ellipsoid model. To reduce any errors arising from the remote sensing data source, Sentinel-1 SAR data was acquired at a time corresponding to the acquisition time of the Sentinel-2 MSI data (Table 1). Sentinel-1 SAR data pre-processing included masking, orbit correction, thermal noise removal, radiometric calibration, terrain correction, and format conversion of the backscattering coefficients, all of which were done with the assistance of SNAP 8.0. Then, a SRTM digital elevation model was utilized with pixel spacing resampled to 10 m by cubic convolution.

In addition to the original remote sensing images, 10 remote sensing indices, 32 texture features, and three topographical features were derived from multi-source data (

Table 2. Candidate variables for image classification and AGB estimation.

Mission	Variable Type		Variable Name	Number
Sentinel-2	Sentinel-2 spectral reflectance		Blue, Green, Red, VRE1, VRE2, VRE3, NIR, N_NIR, SWIR1, SWIR2	10
	Remote sensing index	Vegetation index	RVI, DVI, PVI, NDVI, SAVI, GEMI, S2REP, MTCI	10
		Others	MNDWI, NDBI
	Sentinel-2 textural feature		PCi_MEA, PCi_VAR, PCi_HOM, PCi_CON, PCi_DIS, PCi_ENT, PCi_SEC, PCi_COR	16
Sentinel-1	Sentinel-1 backscattering coefficient		VV, VH	2
	Sentinel-1 texture feature		PP_MEA, PP_VAR, PP_HOM, PP_CON, PP_DIS, PP_ENT, PP_SEC, PP_COR	16
SRTM	Topographical feature		DEM, Slope, Aspect	3

Note: VRE1, VRE2, and VRE3 correspond to the three red-edge bands (~704 nm, ~740 nm, and ~783 nm) of the Sentinel-2 MSI data, and N_NIR represents the near-infrared band (~833 nm). PCi represents band i (1, 2) of the results of principal component analysis of Sentinel-2 MSI data, PP represents the two bands (VV, VH) of Sentinel-1 SAR data; and MEA, VAR, HOM, CON, DIS, ENT, SEC, and COR in order correspond to texture features: mean, variance, homogeneity, contrast, dissimilarity, entropy, angular second moment, and correlation.

). Based on the coefficient of variation (CV) analysis [48] and fineness judgment of the image texture (see Figure A1 and Figure A2, Appendix A), texture feature analysis was performed using a gray level co-occurrence matrix (GLCM) with an optimal window size of 5 × 5 pixels, a shift step of 1 in the X and Y directions, and a gray scale compression selection of 32 levels.

2.3. Urban Green Space Classification

The prime, core part of object-oriented classification is image segmentation, which generates internally “homogeneous and cohesive” polygonal objects that directly impact the following classification process. Image “under-segmentation” or “over-segmentation” can be easily caused by the use of a single segmentation scale, so multi-scale segmentation makes it possible to perform reasonable segmentation of all features [49]. The affiliation function algorithm, which is frequently applied in object-oriented classification, was created based on fuzzy functions and classifies images via the formation of logical functions [50]. Thereby, the distinction between the three major land cover types (urban green spaces, urban impervious surfaces, and water bodies) was first achieved using the above approaches. The strategies for image segmentation were as follows: (1) weights for the red-edge bands (bands 5, 6, and 7) from Sentinel-2 MSI data were taken as 1.20, because they reflect monitored vegetation growth, and the other multispectral bands were weighted as 1; (2) the corresponding optimal scale parameters for urban green spaces, urban impervious surfaces, and water bodies were set to 25, 30, and 50, respectively, due to their different spatial connectivity and heterogeneity; (3) the color parameter was set to 0.80, the shape parameter was set to 0.20, and the compactness and smoothness parameters were both 0.50.

Machine learning methods were used to further classify UGS into three blocks of broadleaf forest, coniferous forest, and shrub-grass vegetation in the object-oriented classification Classifier module of eCognition 9.1. Multiple image features (Table 2) were combined to establish the classification rules. Since the relationships between different vegetation patches in urban ecosystems are highly complex and nonlinear, random forest (RF) and k-nearest neighbor (KNN), both of which can be used for multi-class and nonlinear classification, were the preferred classifiers tested in this study. RF, which combines Bagging integrated learning theory with a random subspace method, is an integrated algorithm, based on decision trees, with high accuracy in the application of massive data-based fast image classification [51,52,53]. KNN is also an effective supervised machine learning method, a lazy learning algorithm that allows the class of each target object to be determined by the k nearest training sample objects in the feature space, according to the Euclidean distance (ED) minimum or majority voting criterion [54]. The commonly used function for weighting distance is:

$${\mathrm{w}}_{\left({\mathrm{T}}_{\mathrm{i}}\right)}{}_{{}_{\mathrm{t}}}=\frac{1}{{\mathrm{d}}_{{{}_{\left({\mathrm{T}}_{\mathrm{i}}\right)}}_{{}_{\mathrm{t}}}}^2}/{\displaystyle \sum _{\mathrm{j}=1}^{\mathrm{k}}\frac{1}{{\mathrm{d}}_{{{}_{\left({\mathrm{T}}_{\mathrm{j}}\right)}}_{{}_{\mathrm{t}}}}^2}}$$

(1)

where ${\mathrm{d}}_{{{}_{\left({\mathrm{T}}_{\mathrm{i}}\right)}}_{{}_{\mathrm{t}}}}$ is the distance between the target object T and the training object T_i in the feature space, and k is the number of closest units.

The confusion matrix was used to evaluate the results of UGS classification with the indexes of producer accuracy (PA), user accuracy (UA), overall classification accuracy (OA), and overall kappa index (KIA).

2.4. Modeling

2.4.1. Correlation Analysis

There were three candidate variable combination scenarios for modeling: (1) variables mainly associated with Sentinel-2 MSI data (n = 37, 10 original bands, eight vegetation indices, 16 texture features, and three topographical features); (2) variables mainly associated with Sentinel-1 SAR data (n = 21, two original bands, 16 texture features, and three topographical features); (3) variables derived from multi-source data (n = 55, 12 original bands, eight vegetation indices, 32 texture features, and three topographical features).

To determine how candidate variables affect AGB of each vegetation type prior to modeling, we performed Pearson’s correlation analysis of these variables with measured biomass data.

2.4.2. Stepwise Regression Modeling

Stepwise regression creates an “optimal” multiple linear regression equation. It sequentially introduces variables into the model and tests the significance of the forward variables and the newly introduced ones on the dependent variable using stepwise calculations, to determine the appropriateness of introducing new independent variables [55]. The F-statistic was used to determine whether the original variables became insignificant after new variables were introduced (the thresholds for entry and removal were set at 0.05 and 0.10, respectively). To improve model significance, attention needs to be paid to the possibility of multi-collinearity among the independent variables [56]. The independent variables with high covariance were excluded by a variance inflation factor (VIF) [57] with the exclusion condition of VIF ≥ 10. We used the stepwise regression model for estimating AGB defined as:

$$\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}}\times {\mathrm{x}}_{\mathrm{i}}+\mathrm{q}}$$

(2)

where $\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}}$ is the simulated value of above-ground biomass, ${\mathrm{p}}_{\mathrm{i}}$ is the regression coefficient of variable i, n is the number of modeling variables, and $\mathrm{q}$ is the regression intercept.

2.4.3. Model Accuracy Assessment

Seventy-five percent of the calculated quadrat biomass values were selected as training data and the remaining were used to estimate model performance. The coefficient of determination (R²) and the root mean square error (RMSE) were applied to assess the credibility of the biomass estimation model.

$${\mathrm{R}}^2={{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}}-\overline{\mathrm{A}\mathrm{G}\mathrm{B}}\right)}}^2/{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{d}}-\overline{\mathrm{A}\mathrm{G}\mathrm{B}}\right)}}^2$$

(3)

$$\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{\frac{1}{\mathrm{n}}{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left(\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{d}}-\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}}\right)}}^2}$$

(4)

where $\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}}$ is the simulated value of above-ground biomass, $\mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{e}\mathrm{d}}$ is the measured value of above-ground biomass, $\overline{\mathrm{A}\mathrm{G}\mathrm{B}}$ is the average of measured values, and n is the number of field data quadrats. Figure 2 is a technical flowchart of the methodology in this study.

3. Results

3.1. Urban Green Space Type Identification

After combining multi-features from remote sensing images, green space coverage maps in urban areas were derived using RF and KNN algorithms and the results are shown in Figure 3. As can be seen in

Table 3. Confusion matrix.

		RF			KNN
User\ Reference Class	Broadleaf Forest	Coniferous Forest	Shrub-Grass Vegetation	Broadleaf Forest	Coniferous Forest	Shrub-Grass Vegetation
Broadleaf forest	604	87	28	587	108	33
Coniferous forest	83	751	0	91	728	0
Shrub-grass vegetation	45	0	240	51	2	208
Unclassified	4	0	0	7	0	27
PA (%)	82.07	89.62	89.55	79.76	86.87	77.61
UA (%)	84.01	90.05	84.21	80.63	88.89	79.69
OA (%)		86.59			82.68
KIA		0.78			0.72

, the overall Kappa coefficients for the UGS classification results by RF (0.78) and KNN (0.72) algorithms both exceeded 0.70, thereby the classification results of UGS were in good agreement with the actual situation. Nevertheless, in terms of overall classification accuracy, the RF (86.59%) results were ~4% greater than KNN (82.68%). Moreover, the confidence level for coniferous forest (producer accuracy of 89.62%) was higher than that of shrub-grass vegetation (producer accuracy of 89.55%) and broadleaf forest (producer accuracy of 82.07%) in the RF classification results.

Finally, object-oriented classification of green space in the study area during the 2016, 2018, and 2021 growing seasons was carried out by training the RF algorithm (Figure 4). Some details of the post-classification results, versus the ground-truth situation, are displayed in Figure 5. Greater insight into the UGS composition structure and change trend is provided by an area statistical table (

Table 4. Areas for different types of UGS in the study area.

Year	Area (km2)				Percentage of the Study Area (%)
	Broadleaf Forest	Coniferous Forest	Shrub-Grass Vegetation	Total
2016	68.84	24.99	159.28	253.10	48.01%
2018	56.28	22.97	192.63	271.88	51.57%
2021	85.38	17.41	142.75	245.54	46.57%

). Shrub-grass vegetation (~58.14% for 2021) was most widespread in approximately half of the UGS in the study area, followed by broadleaf forests (~34.77% for 2021) and finally coniferous forests (~7.09% for 2021). During the study period, the area of broadleaf forest decreased and then increased, the area of coniferous forest decreased slightly, and the area of shrub-grass vegetation first increased and then decreased in line with the general trend of UGS change.

3.2. Above-Ground Biomass Modeling

3.2.1. Correlation Analysis

Twenty-six, twenty-three, and eighteen variables were significantly correlated with the measured AGB in the broadleaf, coniferous, and shrub-grass vegetation quadrats, respectively (Figure 6). For broadleaf forests, 13 variables were positively correlated with biomass, of which the highest positive correlation was VH (0.48), whereas most vegetation indices were negatively correlated with biomass (−0.44 to −0.37) and the highest negative correlation was SWIR1 (−0.47). For coniferous forests, 10 variables were positively correlated with biomass, and the SAVI and GEMI vegetation indices were positively correlated with biomass. For shrub-grass vegetation, VH (0.50), DEM (0.44), and Slope (0.44) were most strongly positively correlated, and SWIR1 (−0.41), Red (−0.44), and SWIR2 (−0.51) were the most strongly negatively correlated.

3.2.2. Stepwise Regression Models

For different green space vegetation types, AGB models, based on Sentinel-2 MSI data, Sentinel-1 SAR data, and multi-source remote sensing data were developed accordingly. The type-specific vegetation AGB estimation models are shown in

Table 5. Stepwise regression models for AGB of different types of UGS.

Dataset	Type	Model	Adj-Rjm2
Sentinel-2	BF	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{s}\mathrm{i}\mathrm{m}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}}& =-2274.36\times \mathrm{S}\mathrm{A}\mathrm{V}\mathrm{I}+1.07\times \mathrm{V}\mathrm{R}\mathrm{E}3-781.53\hfill \\ & \times \mathrm{P}\mathrm{C}1\_\mathrm{M}\mathrm{E}\mathrm{A}-1708.51\times \mathrm{P}\mathrm{C}1\_\mathrm{D}\mathrm{I}\mathrm{S}-624.45\hfill \\ & \times \mathrm{P}\mathrm{C}1\_\mathrm{C}\mathrm{O}\mathrm{R}-1403.07\times \mathrm{P}\mathrm{C}2\_\mathrm{H}\mathrm{O}\mathrm{M}+240.33\hfill \\ & \times \mathrm{P}\mathrm{C}2\_\mathrm{C}\mathrm{O}\mathrm{R}+3490.70\hfill \end{array}$	0.37
	CF	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =-12.48\times \mathrm{Slope}+5.26\times \mathrm{GEM}I+1.74\times \mathrm{Red}\hfill \\ & -9024.21\times \mathrm{PC}1\_VAR-2828.51\times \mathrm{PC}1\_\mathrm{SEC}\hfill \\ & -605.03\times \mathrm{PC}2\_\mathrm{ENT}-886.32\times \mathrm{PC}2\_\mathrm{COR}\hfill \\ & +3295.78\hfill \end{array}$	0.58
	SGV	$\begin{array}{cc}\hfill \mathrm{AG}{\mathrm{B}}_{\mathrm{simulated}}& =2.43\times \mathrm{S}2\mathrm{R}\mathrm{E}\mathrm{P}+1.41\times \mathrm{D}\mathrm{E}\mathrm{M}+5.89\times \mathrm{S}\mathrm{l}\mathrm{o}\mathrm{p}\mathrm{e}\hfill \\ & -57.89\times \mathrm{P}\mathrm{C}2\_\mathrm{C}\mathrm{O}\mathrm{N}-1751.94\times \mathrm{P}\mathrm{C}1\_\mathrm{H}\mathrm{O}\mathrm{M}\hfill \\ & -237.54\times \mathrm{P}\mathrm{C}1\_\mathrm{E}\mathrm{N}\mathrm{T}-22.42\hfill \end{array}$	0.61
Sentinel-1	BF	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =92.17\times \mathrm{V}\mathrm{H}-95.13\times \mathrm{V}\mathrm{V}+4729.59\hfill \\ & \times \mathrm{V}\mathrm{H}\_\mathrm{H}\mathrm{O}\mathrm{M}-2117.72\times \mathrm{V}\mathrm{H}\_\mathrm{S}\mathrm{E}\mathrm{C}+284.52\hfill \\ & \times \mathrm{V}\mathrm{V}\_\mathrm{M}\mathrm{E}\mathrm{A}-12.20\times \mathrm{D}\mathrm{E}\mathrm{M}-6277.24\hfill \end{array}$	0.54
	CF	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =0.92\times \mathrm{Aspect}-313.11\times \mathrm{V}\mathrm{V}\_\mathrm{C}\mathrm{O}\mathrm{R}\hfill \\ & -193.94\times \mathrm{V}\mathrm{H}\_\mathrm{C}\mathrm{O}\mathrm{N}+2551.59\times \mathrm{V}\mathrm{V}\_\mathrm{S}\mathrm{E}\mathrm{C}\hfill \\ & +1239.67\times \mathrm{V}\mathrm{V}\_\mathrm{E}\mathrm{N}\mathrm{T}-1896.82\hfill \end{array}$	0.55
	SGV	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =5.37\times \mathrm{V}\mathrm{H}+5.82\times \mathrm{V}\mathrm{V}+2.65\times \mathrm{Slope}\hfill \\ & +65.90\times \mathrm{V}\mathrm{H}\_\mathrm{E}\mathrm{N}\mathrm{T}-28.19\times \mathrm{V}\mathrm{V}\_\mathrm{M}\mathrm{E}\mathrm{A}\hfill \\ & -80.21\times \mathrm{V}\mathrm{H}\_\mathrm{C}\mathrm{O}\mathrm{N}+580.33\hfill \end{array}$	0.78
Sentinel-2 and Sentinel-1	BF	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =-0.08\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}1+23.27\times \mathrm{V}\mathrm{H}-7.97\times \mathrm{S}2\mathrm{R}\mathrm{E}\mathrm{P}\hfill \\ & +25.85\times \mathrm{M}\mathrm{T}\mathrm{C}\mathrm{I}-36.70\times \mathrm{V}\mathrm{H}\_\mathrm{M}\mathrm{E}\mathrm{A}\hfill \\ & +86.19\times \mathrm{P}\mathrm{C}1\_\mathrm{V}\mathrm{A}\mathrm{R}-57.67\times \mathrm{P}\mathrm{C}2\_\mathrm{V}\mathrm{A}\mathrm{R}\hfill \\ & -138.59\times \mathrm{P}\mathrm{C}2\_\mathrm{S}\mathrm{E}\mathrm{C}+50.66\times \mathrm{P}\mathrm{C}2\_\mathrm{C}\mathrm{O}\mathrm{R}\hfill \\ & +21.76\times \mathrm{V}\mathrm{V}\_\mathrm{M}\mathrm{E}\mathrm{A}+6554.20\hfill \end{array}$	0.78
	CF	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =-13.58\times \mathrm{Slope}-2.70\times \mathrm{Blue}+1018.04\hfill \\ & \times \mathrm{V}\mathrm{H}\_\mathrm{S}\mathrm{E}\mathrm{C}+0.53\times \mathrm{Aspect}+203.72\hfill \\ & \times \mathrm{P}\mathrm{C}2\_\mathrm{C}\mathrm{O}\mathrm{R}+0.57\times \mathrm{S}\mathrm{W}\mathrm{I}\mathrm{R}2-610.44\hfill \\ & \times \mathrm{V}\mathrm{H}\_\mathrm{C}\mathrm{O}\mathrm{N}+862.00\times \mathrm{V}\mathrm{V}\_\mathrm{S}\mathrm{E}\mathrm{C}+620.30\hfill \\ & \times \mathrm{V}\mathrm{H}\_\mathrm{E}\mathrm{N}\mathrm{T}+961.30\times \mathrm{V}\mathrm{V}\_\mathrm{C}\mathrm{O}\mathrm{N}+1280.68\hfill \end{array}$	0.69
	SGV	$\begin{array}{cc}\hfill \mathrm{A}\mathrm{G}{\mathrm{B}}_{\mathrm{simulated}}& =7.68\times \mathrm{V}\mathrm{H}-74.91\times \mathrm{V}\mathrm{H}\_\mathrm{C}\mathrm{O}\mathrm{N}+\hfill \\ & \times \mathrm{V}\mathrm{H}\_\mathrm{E}\mathrm{N}\mathrm{T}+1.32\times \mathrm{D}\mathrm{E}\mathrm{M}-5.60\times \mathrm{M}\mathrm{T}\mathrm{C}\mathrm{I}\hfill \\ & -12.02\times \mathrm{V}\mathrm{V}\_\mathrm{M}\mathrm{E}\mathrm{A}+275.82\hfill \end{array}$	0.80

Note: BF, CF, and SGV, refer to broadleaf forest, coniferous forest, and shrub-grass vegetation.

.

3.2.3. Accuracy Assessment

The accuracy evaluation results from the AGB estimation models are shown in Figure 7. The R² values of all models exceeded 0.50. Among the biomass models, the poorest performance was found in the broadleaf forest biomass estimation model, based on Sentinel-2 MSI data, with both the lowest R² value (0.54) and the highest RMSE (54.88 tons/ha). The coniferous forest biomass estimation model, with a combination of Sentinel-2 and Sentinel-1 data, had the highest R² value (0.77). The combination of Sentinel-2 and Sentinel-1 data used to estimate the biomass of shrub-grass vegetation had the lowest RMSE (15.05 tons/ha).

3.3. Above-Ground Biomass Estimation

The optimal models built by combining the use of Sentinel-2 MSI data and Sentinel-1 SAR data were selected from the model estimation results to estimate AGB of UGS in the study area at different periods. It was derived by overlaying the AGB maps for the three green space types. Figure 8 shows that AGB for Xuzhou UGS decreased generally from southwest to northeast, with a distribution range of 0–160 tons/ha.

The total AGB values in 2016, 2018, and 2021 were estimated to be 1.40 × 10⁶ tons, 1.47 × 10⁶ tons, and 1.61 × 10⁶ tons, respectively. The AGB for broadleaf forests in 2016, 2018, and 2021 accounted for 57.76%, 56.75%, and 70.86% of the total in the corresponding periods. More importantly, the average value for UGS in 2021 increased by about 10 tons/ha compared to 2016 (Figure 9).

4. Discussion

The combination of object-oriented classification and machine learning methods using high-resolution multispectral or hyperspectral image-based classification proved to distinguish urban forest from other urban landscapes [58,59,60]. To further distinguish diverse urban vegetation components, such as different species of trees and shrub-grass vegetation, in this study, we fused optical and SAR images for green space type identification and further confirmed that urban green space (UGS) classification based on image objects using an RF algorithm was more effective than the typical KNN. This method had the highest classification accuracy for coniferous forests (producer accuracy of 89.62%), which may have been due to concentrated and continuous distribution and orderly spatial arrangement of coniferous forests with unique stand structure characteristics.

The modeling pre-selection variable screening showed that for all three types of green spaces, the spectral reflectance of Sentinel-2 MSI data was mostly negatively correlated with above-ground biomass (AGB) of green space. In addition, vegetation indices and texture features influenced AGB estimation, which has been reported in the results of other studies [61,62]. A strong correlation was found between the variable aspect (0.54) and the AGB of coniferous forests. This was because coniferous forests in Xuzhou are typically located in mountainous areas and the specific geographical location influences the length of sunlight available to trees, soil fertility, etc. Thus, it should be particularly noted that topographic features were likewise key influencing variables in urban areas with undulating landscapes.

Previous studies have shown that vegetation AGB can be effectively and robustly estimated using optical or radar remote sensing images [63,64,65,66]. Our study equally found that AGB of coniferous forest could be estimated using Sentinel-2 MSI data (R² = 0.69, RMSE = 22.83 tons/ha). Moreover, urban shrub-grass vegetation biomass could be better estimated using Sentinel-1 (R² = 0.69, RMSE = 15.72 tons/ha) than Sentinel-2 data (R² = 0.55, RMSE = 20.30 tons/ha) in this study, in contrast with the findings on Mediterranean shrubland biomass by Chang et al. (R² < 0.60 for Sentinel-1 model, but R² = 0.72 for Sentinel-2 model) [67]. Besides this, compared with a quantitative study of Xuzhou broadleaf forest biomass estimation using only Sentinel-2 MSI data by Li et al. (R² = 0.73, RMSE = 45.56 tons/ha) [62], we proved that the additional use of Sentinel-1 SAR data showed better accuracy (R² = 0.75, RMSE = 35.95 tons/ha). More importantly, using both data sources provided better estimates of AGB compared to modeling with Sentinel-2 MSI data or Sentinel-1 SAR data alone. The results of a study by Nuthammachot et al. also revealed higher accuracy of the model (R² = 0.79, RMSE = 22.98 tons/ha) when combining Sentinel-1 and Sentinel-2 data than modeling with Sentinel-2 data only (R² = 0.77, RMSE = 24.19 tons/ha) for predicting forest biomass [68]. Meanwhile, Wang et al. demonstrated the potential of the simultaneous use of HJ-1B and RadarSat-2 data to estimate AGB in the desert grasslands of Ningxia, China (R² = 0.71, RMSE = 14.20 kg/hm²), which outperformed the AGB estimation based on the NDVI index of HJ-1B (R² = 0.27, RMSE = 20.58 kg/hm²) [69]. Therefore, we argue that combining the strengths of optical remote sensing data and radar data enables better quantitative estimation of AGB of UGS, an inference that agrees with several studies [70,71,72,73].

Despite the overall decline in the total area of UGS in Xuzhou, the AGB showed an opposite increasing trend during the study period, mainly due to the increase in the area of broadleaf forests with high biomass attributes. Moreover, the average AGB of broadleaf forests (133.71 tons/ha in 2021) surpassed that of coniferous forests (92.13 tons/ha in 2021), as evidenced by the findings of Xing et al. [74]. There was no clear annual variability in the average AGB of shrub-grass vegetation, with a mean value as low as 21.81 tons/ha from 2016 to 2021; nevertheless, its contribution to AGB in green areas should not be underestimated, as its coverage area was considerably higher than the other two types of vegetation. Therefore, for the future construction of UGS systems, adequate broadleaf forests ought to be reserved, while giving high value to the coordination arrangement of broadleaf forests, coniferous forests, and shrub-grass vegetation. This will improve green space maintenance and overall environmental outcomes.

At the same time, there are still some limitations. The classification results for broadleaf forests had relatively low, but acceptable, accuracy (producer accuracy of 82.07%), which was caused by the misclassification of objects belonging to broadleaf forests as coniferous forests. Consequently, it remains to be tested whether the additional use of vegetation phenological characteristics and thermal properties from remote sensing images would help reduce misclassification. This study also verified that the backscatter coefficient and texture features of Sentinel-1 alone can hardly reflect the full variability of AGB in green spaces. Therefore, it would be more interesting if richer image features were derived to further test the capabilities of Sentinel-1 SAR data for AGB estimation.

5. Conclusions

This study demonstrated the feasibility of applying a combination of Sentinel-2 MSI data and Sentinel-1 SAR data to map urban green space and estimate above-ground biomass. The main findings of the study are as follows:

Based on Sentinel satellite series high-resolution remote sensing images, the use of object-oriented classification and random forest algorithm is a reliable method to identify urban green space types.

The integration of Sentinel-2 MSI data and Sentinel-1 SAR data improved the accuracy of a stepwise regression model for above-ground biomass in urban green space. Model prediction accuracy when using combined multi-source data was, from highest to lowest, coniferous forest (R² = 0.77, RMSE = 18.39 tons/ha), shrub-grass vegetation (R² = 0.76, RMSE = 15.05 tons/ha), and broadleaf forest (R² = 0.75, RMSE = 35.95 tons/ha).

From 2016 to 2021, the above-ground biomass of urban green space in Xuzhou increased, most prominently due to broadleaf forests. The highest average above-ground biomass was derived from broadleaf forests and the lowest from shrub-grass vegetation.

Altogether, the combination high-resolution Sentinel-2 MSI data and Sentinel-1 SAR data provides a practical solution for estimating above-ground biomass of urban green space in other regions. Much effort has been made in this study to consider the urban green space dynamics in Xuzhou in the last five years. Further research on natural environmental factors, such as soil physical and chemical properties, is needed to improve the above-ground biomass distribution and mapping accuracy of urban green space in Xuzhou.