Spatial Pattern Analysis of the Ecosystem Services in the Guangdong-Hong Kong-Macao Greater Bay Area Using Sentinel-1 and Sentinel-2 Imagery Based on Deep Learning Method

: Assessment of ecosystem services supply, demand, and budgets can help to achieve sustainable urban development. The Guangdong-Hong Kong-Macao Greater Bay Area, as one of the most developed megacities in China, sets up a goal of high-quality development while fostering ecosystem services. Therefore, assessing the ecosystem services in this study area is very important to guide further development. However, the spatial pattern of ecosystem services, especially at local scales, is not well understood. Using the available 2017 land cover product, Sentinel-1 SAR and Sentinel-2 optical images, a deep learning land cover mapping framework integrating deep change vector analysis and the ResUnet model was proposed. Based on the produced 10 m land cover map for the year 2020, recent spatial patterns of the ecosystem services at different scales (i.e., the GBA, 11 cities, urban–rural gradient, and pixel) were analyzed. The results showed that: (1) Forest was the primary land cover in Guangzhou, Huizhou, Shenzhen, Zhuhai, Jiangmen, Zhaoqing, and Hong Kong, and an impervious surface was the main land cover in the other four cities. (2) Although ecosystem services in the GBA were sufﬁcient to meet their demand, there was undersupply for all the three general services in Macao and for the provision services in Zhongshan, Dongguan, Shenzhen, and Foshan. (3) Along the urban–rural gradient in the GBA, supply and demand capacity showed an increasing and decreasing trend, respectively. As for the city-level analysis, Huizhou and Zhuhai showed a ﬂuctuation pattern while Jiangmen, Zhaoqing, and Hong Kong presented a decreasing pattern along the gradient. (4) Inclusion of neighborhood landscape led to increased demand scores in a small proportion of impervious areas and oversupply for a very large percent of bare land.


Introduction
Ecosystem services, defined as benefits that humans can obtain from the ecosystem, have been increasingly investigated and considered in both research and policy communities due to raising awareness among decision-makers and the public [1][2][3]. The maintenance of ecosystem services for enhancing human well-being was included in a set of goals of the 2050 Vision of the Strategic Plan for Biodiversity of the Convention on Biological Diversity [4]. In addition, as summarized in [5], China has established several relevant policies, for example, a national Forest Eco-Compensation Fund system [6] and ecological redline policy [7]. Simultaneously, due to urbanization and population growth, ecosystem services are continuously in decline. As indicated in [8], the global loss of ecosystem services due to land use changes from 1997 to 2011 was estimated at $US 4.3-20.2 trillion/year. Therefore, the assessment of ecosystem services is very important for decision-makers to achieve sustainable urban development.
The ecosystem services can be quantified in semi-quantitatively (e.g., capacity scores) [9,10], monetary units (e.g., yuan and dollar) [11,12], and biophysical units (e.g., biomass) [13,14]. Among the different available methods, expert-based ecosystem services scores are considered to be flexible and adaptable [1]. One such method is the "capacity matrix approach" linking supply, demand, and budget capacity to land use/land cover types [15], which was referred to as "one of the most popular ecosystem services assessment techniques today" in [16]. It has been widely applied in many studies. For example, ecosystem services capacity scores for each land use type and their total supply, demand, and budgets in Shenzhen were mapped and monitored in [9]. In [17], the time-varying characteristics of different ecosystem services of China on three scales: the total quantity change on a national scale, the time-varying trends on a provincial scale, and the change rates on a city scale were analyzed. An assessment of the spatial pattern of ecosystem services across Europe on a 1 km 2 grid was presented [18]. Although extensive studies at different spatial scales ranging from regional [9] to national [17], and even to continental [18,19] scales have been conducted, the spatial pattern of ecosystem services at higher spatial resolution (10 m or less) is rarely concerned. Unfortunately, relevant information for local-scale decision-making is still lacking.
Since remote sensing can provide data at a low cost, it is appropriate for ecosystem services assessment at higher spatial resolution. The free availability of medium-high spatial resolution imagery, such as Sentinel-1 SAR and Sentinel-2 optical data, can facilitate the mapping of finer scale (10 m) ecosystem services. Compared with coarse or medium spatial resolution data, they have the advantage of improved spatial details for the identification of micro-scale land cover/land use (e.g., small vegetation patch in urban area). In addition, the fusion of Sentinel-1 and Sentinel-2 has been proved to be effective in improving land use/land cover mapping accuracies [20,21] since complemented backscatter and spectral information are utilized [22]. In this way, the reliability of ecosystem services analysis is guaranteed. For the land use/land cover mapping techniques, advanced artificial neural networks, especially deep learning models, have gained increased attention in the remote sensing domain owing to the end-to-end nature (i.e., learn to discriminate features designed for image classification and associated classifier jointly [23]). One of the obstacles affecting the performance of deep learning networks is the limited training samples problem [24], while training samples collection work is expensive and time-consuming. Therefore, a deep learning-based land use/land cover mapping framework that can fully integrate Sentinel SAR-optical information and handle the training samples lacking problem should be established.
Based on the two aforementioned aspects, the overarching goal of this paper is to develop a land cover mapping architecture using Sentinel-1 and Sentinel-2 imagery based on the deep learning method and investigate the spatial pattern of ecosystem services at a finer spatial scale. To this end, we take the Guangdong-Hong Kong-Macao Greater Bay Area (GBA), one of the most important economic regions in China, as our study area. Firstly, a deep learning method combing deep change vector analysis and model fine-tuning was proposed for land cover mapping. Then, the ecosystem service values were quantified by linking land cover types to supply and demand matrixes. Finally, ecosystem services scores at different spatial scales (i.e., the GBA, different cities, buffer zones, and pixels) were investigated and analyzed.

Study Areas
The GBA, situated at the Pearl River Delta, is composed of 9 mainland cities and 2 special administrative regions, including Guangzhou, Shenzhen, Zhuhai, Foshan, Huizhou, Dongguan, Zhongshan, Jiangmen, Zhaoqing, Hong Kong, and Macao ( Figure 1). With 56,000 km 2 coverage (about 0.6% of China's land area), the GBA created a gross domestic product (GDP) of more than $11 trillion in 2019, contributing about 1/7 of the national GDP. In February 2019, the Chinese government released the Outline Development Plan for the Guangdong-Hong Kong-Macao Greater Bay Area. It clarified the strategic positioning of a vibrant world-class agglomeration, an international technology and innovation center, an important support pillar for the Belt and Road Initiative, a showcase for in-depth cooperation between the mainland, Hong Kong, and Macao, and a high-quality living circle [25]. Under the collaborative development of the economy and the environment, balancing urban development and ecological conservation will become increasingly prominent. Therefore, analysis of ecosystem service value in the GBA is fundamental to guide sustainable urban development. With 56,000 km 2 coverage (about 0.6% of China's land area), the GBA created a gross domestic product (GDP) of more than $11 trillion in 2019, contributing about 1/7 of the national GDP. In February 2019, the Chinese government released the Outline Development Plan for the Guangdong-Hong Kong-Macao Greater Bay Area. It clarified the strategic positioning of a vibrant world-class agglomeration, an international technology and innovation center, an important support pillar for the Belt and Road Initiative, a showcase for in-depth cooperation between the mainland, Hong Kong, and Macao, and a high-quality living circle [25]. Under the collaborative development of the economy and the environment, balancing urban development and ecological conservation will become increasingly prominent. Therefore, analysis of ecosystem service value in the GBA is fundamental to guide sustainable urban development.

Datasets
In this study, the 10 m land cover product FROM-GLC10 [26] and Sentinel-1 and Sentinel-2 remote sensing datasets were used. FROM-GLC10 is a 10 m global land cover product of the year 2017 and has 10 land cover types, including cropland, forest, grassland, shrubland, wetland, water, tundra, impervious surface, bare land, and snow/ice. Tundra and snow/ice are not included in our study area. Wetland is merged into grassland since they show very similar spectral signature, and some classified wetland pixels are found to be grassland. Therefore, 7 land cover types (i.e., cropland, forest, grassland, shrubland, water, impervious surface, and bare land) were employed as our classification system. The remote sensing datasets composed of Sentinel-1 SAR ground range detected (GRD) images and Sentinel-2 optical images were obtained from Sentinels Scientific Data Hub (https://scihub.copernicus.eu/; accessed on 18 February 2021). The pre-processing steps of Sentinel-1 GRD data, including calibration, speckle-filter, and terrain correction, were implemented by the European Space Agency's (ESA) Sentinel Application Platform (SNAP), generating geocoded VV and VH backscattering coefficients with the spatial resolution of 10 m. Sentinel-2 satellites can provide four 10 m bands, six 20 m bands, and three 60 m bands. Sentinel-2 datasets were processed into atmospherically corrected products by the ESA's Sen2Cor algorithm [27]. Only four 10 m bands (i.e., blue, green, red, and near-infrared bands) of Sentinel-2 data were used. Multi-temporal Sentinel optical-SAR

Datasets
In this study, the 10 m land cover product FROM-GLC10 [26] and Sentinel-1 and Sentinel-2 remote sensing datasets were used. FROM-GLC10 is a 10 m global land cover product of the year 2017 and has 10 land cover types, including cropland, forest, grassland, shrubland, wetland, water, tundra, impervious surface, bare land, and snow/ice. Tundra and snow/ice are not included in our study area. Wetland is merged into grassland since they show very similar spectral signature, and some classified wetland pixels are found to be grassland. Therefore, 7 land cover types (i.e., cropland, forest, grassland, shrubland, water, impervious surface, and bare land) were employed as our classification system. The remote sensing datasets composed of Sentinel-1 SAR ground range detected (GRD) images and Sentinel-2 optical images were obtained from Sentinels Scientific Data Hub (https://scihub.copernicus.eu/; accessed on 18 February 2021). The pre-processing steps of Sentinel-1 GRD data, including calibration, speckle-filter, and terrain correction, were implemented by the European Space Agency's (ESA) Sentinel Application Platform (SNAP), generating geocoded VV and VH backscattering coefficients with the spatial resolution of 10 m. Sentinel-2 satellites can provide four 10 m bands, six 20 m bands, and three 60 m bands. Sentinel-2 datasets were processed into atmospherically corrected products by the ESA's Sen2Cor algorithm [27]. Only four 10 m bands (i.e., blue, green, red, and near-infrared bands) of Sentinel-2 data were used. Multi-temporal Sentinel optical-SAR images were acquired almost simultaneously around 26 October 2017 and 29 October 2020 to cover the entire study area, referred to as the 'GBA 2017' and 'GBA 2020' images.

Land Cover Mapping
The methodological framework of the deep learning-based land cover mapping is presented in Figure 2, which has three major components: (1) Residual Unet (ResUnet) training with pre-temporal data, (2) sample generation with deep change vector analysis (CVA), and (3) fine-tuning ResUnet with post-temporal data. ResUnet training with pretemporal data refers to model training with 2017 labeled data, serving as feature extraction for the following CVA. Samples generation with deep CVA is to produce training samples for the classification of the GBA 2020 images. Moreover, fine-tuning ResUnet with posttemporal data is to fine-tune the network pre-trained on the GBA 2017 images, making it more adaptable to the 2020 image classification tasks.
deep feature extraction in CVA and image classification. The architecture of the ResUne is shown in Figure 3. ResUnet is an integrated architecture incorporating Unet and dee residual learning, which has three parts: encoding, decoding, and bridge [28]. On the on hand, Unet, a pixel-wise semantic segmentation structure, can combine low-level deta information and high-level semantic information [30], thus enabling promising perfo mance in feature extraction and classification [31]. On the other hand, residual blocks ca ease deep network training by preserving input features through identity skip connec tions [28]. Therefore, ResUnet was used in this study, considering that it is appropriate fo accurate pixel-wise classification and fine-scale spatial pattern analysis. Furthermore, t ensure that information of small targets in the 10 m spatial resolution images can be bette preserved after convolution and pooling layers, while higher abstract features can be cap tured, a 7-level architecture including the total of 15 convolution layers was utilized. A shown in Figure 3, image patches with the size of 512 × 512 × 6 (i.e., four spectral and tw backscatter intensity bands) are fed to the network, and the 512 × 512 × 7 classificatio probability map is produced at the output layer.  Sentinel-1 optical/Sentinel-2 SAR images with the spatial resolution of 10 m are rich in spectral/backscatter and spatial information. Inspired by the effective capturing of this information by ResUnet in optical [28] and SAR [29] images, it was used for effective feature extraction and following classification. In this study, ResUnet was adopted for both deep feature extraction in CVA and image classification. The architecture of the ResUnet is shown in Figure 3. ResUnet is an integrated architecture incorporating Unet and deep residual learning, which has three parts: encoding, decoding, and bridge [28]. On the one hand, Unet, a pixel-wise semantic segmentation structure, can combine low-level detail information and high-level semantic information [30], thus enabling promising performance in feature extraction and classification [31]. On the other hand, residual blocks can ease deep network training by preserving input features through identity skip connections [28]. Therefore, ResUnet was used in this study, considering that it is appropriate for accurate pixel-wise classification and fine-scale spatial pattern analysis. Furthermore, to ensure that information of small targets in the 10 m spatial resolution images can be better preserved after convolution and pooling layers, while higher abstract features can be captured, a 7-level architecture including the total of 15 convolution layers was utilized. As shown in Figure 3, image patches with the size of 512 × 512 × 6 (i.e., four spectral and two backscatter intensity bands) are fed to the network, and the 512 × 512 × 7 classification probability map is produced at the output layer.
For sample generation, unchanged locations between 2017 and 2020 were identified through deep CVA [32], and the corresponding unchanged 2020 images were labeled according to land cover types in the FROM-GLC10 map. Deep CVA starts from the extraction of deep features that can model low-level spatial context and high-level semantic information. Deep CVA proposed by Saha et al. [31] employed a pre-trained CNN model trained on aerial RGB image bands for deep feature extraction. Since it was trained on different data sets and tasks, extracted deep features are not optimal. Therefore, in our study, firstly, the ResUnet model was trained with stacked Sentinel 1 and 2 image bands in 2017 as inputs and FROM-GLC10 map as labels. With the trained ResUnet model, the multi-temporal deep features were obtained by feeding the GBA 2017 and 2020 images separately as inputs. Subsequently, deep change vectors were generated by subtracting the multi-temporal deep features. Locations with a magnitude of deep change vectors smaller than 10 were identified as unchanged. The threshold value of 10 was determined based on visual inspection of CVA magnitude results. It should be noted that a deliberately selected strict threshold was used to ensure that detected unchanged regions were highly reliable. Finally, with the generated 2020 training samples, the trained ResUnet model can be fine-tuned to suit the classification of the GBA in the year 2020. For sample generation, unchanged locations between 2017 and 2020 were iden through deep CVA [32], and the corresponding unchanged 2020 images were label cording to land cover types in the FROM-GLC10 map. Deep CVA starts from the e tion of deep features that can model low-level spatial context and high-level seman formation. Deep CVA proposed by Saha et al. [31] employed a pre-trained CNN m trained on aerial RGB image bands for deep feature extraction. Since it was train different data sets and tasks, extracted deep features are not optimal. Therefore, i study, firstly, the ResUnet model was trained with stacked Sentinel 1 and 2 image in 2017 as inputs and FROM-GLC10 map as labels. With the trained ResUnet mod multi-temporal deep features were obtained by feeding the GBA 2017 and 2020 im separately as inputs. Subsequently, deep change vectors were generated by subtr the multi-temporal deep features. Locations with a magnitude of deep change v smaller than 10 were identified as unchanged. The threshold value of 10 was determ based on visual inspection of CVA magnitude results. It should be noted that a de ately selected strict threshold was used to ensure that detected unchanged regions highly reliable. Finally, with the generated 2020 training samples, the trained Re model can be fine-tuned to suit the classification of the GBA in the year 2020.

Spatial Pattern Analysis of Ecosystem Services
Assessment of ecosystem services was conducted based on the expert-based s and demand matrixes proposed in [15]. The supply and demand capacity attribu each class are expert evaluation values. In the matrix, capacities of different land

Spatial Pattern Analysis of Ecosystem Services
Assessment of ecosystem services was conducted based on the expert-based supply and demand matrixes proposed in [15]. The supply and demand capacity attributed to each class are expert evaluation values. In the matrix, capacities of different land cover types to supply/demand 22 ecosystem services (three general categories including regulating, provisioning, and cultural services) were provided. A relative scale ranging from 0 to 5 was used to indicate no relevant capacity, low relevant capacity, relevant capacity, medium relevant capacity, high relevant capacity, and very high relevant capacity of supplying/demanding a certain ecosystem service within the particular land cover type. There are 11, 9, and 2 sub-services in regulating, provisioning, and cultural services, respectively. Due to the clarity of the results, only the three general ecosystem categories were considered. With this regard, supply and demand capacity scores were calculated by multiplying the score of each individual ecosystem service for a certain land cover type in the matrix by the area of this land cover type in the region of interest and then summed, and subsequently demand capacity scores were subtracted from the supply capacity scores to derive budget scores ( Table 1). As can be seen in Table 1, the main human-dominated land cover type (i.e., impervious surface) has a higher demand capacity, and the natural land cover types including forest, grassland, shrubland, and water are characterized by higher supply capacity and lower demand capacity. By linking the land cover types in the classification map to corresponding capacity scores, fine-scale ecosystem service supply, demand, and budget maps were generated. Three aspects were considered to analyze the spatial pattern of ecosystem services in the GBA: total ecosystem services values, ecosystem services values along the urban-rural gradient, and ecosystem services balances at the local scale. The total ecosystem services values were investigated for the GBA and the 11 cities. For the ecosystem services values along the urban-rural gradient, buffer zones emanating outward from the urban center (i.e., the centroid of impervious surface area) were constructed. The radiation distances between each buffer zone for the GBA and the 11 cities were set as 10 km and 5 km, respectively. The constructed buffer zones were supposed to cover the GBA and each city. In particular, the buffer zones emanating outward from cities' urban centers were constructed to cover each city, and zone areas that overlapped with neighboring cities were kept. In this way, spatial neighboring effects can be considered and analyzed. The mean values of ecosystem services in each buffer zone were calculated to ensure the evaluation metric independent of the buffer zone size and the shape of cities. Focusing on distance from urban center versus mean values of ecosystem services, their spatial patterns along the urban-rural gradient can be presented. Finally, ecosystem services balances between the demand capacity of a certain pixel and the supply capacity of its neighboring pixels were considered, which was designed to measure the influence of neighborhood landscape on ecosystem services supply and demand balances at a fine spatial scale. In specific, a local neighborhood with a certain circle radius R was defined. Based on the underlying assumption that neighboring pixels that are closer to the center one can provide greater supply capacity, the inverse distance weights were used. On the other hand, since each pixel can provide supply capacity for pixels within the neighborhood, the inverse distance was multiplied by a constant coefficient K to ensure the sum of weight values equal to 1. Ecosystem services balance for each pixel was calculated as follows: where pixels and center denote all pixels and center one in a local neighborhood, Dist i is the distance between pixel i in a local neighborhood and a center one, S and D are the supply and demand capacity scores, respectively. In this study, R was set as 200 m, 500 m, and 1000 m to exploit the effect of multi-scale landscape on ecosystem services balance.

Land Cover Mapping Results
The accuracy of the land cover mapping network around 90% was achieved, which indicated that the proposed framework was able to generate reliable classification results. The land cover mapping result is presented in Figure 4. The impervious surface located in the middle of the GBA and the impervious surface of some cities were spatially connected, for example, Foshan and Guangzhou. Forestland was mainly distributed in the upper-left and upper-right regions of the GBA (i.e., Zhaoqing and Huizhou). Cropland, taking up a small proportion, was mainly distributed in the rural area of Zhaoqing, Huizhou, and Jiangmen. The GBA was rich in natural water resources. A dense river network was staggered in the impervious surface area, and ponds and lakes were intertwined with rivers or distributed in the forest area.

Land Cover Mapping Results
The accuracy of the land cover mapping network around 90% was achieved, which indicated that the proposed framework was able to generate reliable classification results. The land cover mapping result is presented in Figure 4. The impervious surface located in the middle of the GBA and the impervious surface of some cities were spatially connected, for example, Foshan and Guangzhou. Forestland was mainly distributed in the upper-left and upper-right regions of the GBA (i.e., Zhaoqing and Huizhou). Cropland, taking up a small proportion, was mainly distributed in the rural area of Zhaoqing, Huizhou, and Jiangmen. The GBA was rich in natural water resources. A dense river network was staggered in the impervious surface area, and ponds and lakes were intertwined with rivers or distributed in the forest area. The percentages of different land cover types in the GBA and different cities are presented in Figure 5 to reveal land cover composition. Results indicated that forest, cropland, and impervious surface were the three major land cover types in the GBA. As for the eleven cities, forests were the primary land cover in Guangzhou, Huizhou, Shenzhen, Zhuhai, Jiangmen, Zhaoqing, and Hong Kong; and the other four cities (i.e., Dongguan, Foshan, Zhongshan, and Macao) had the largest proportion in impervious surface. Land urbanization that can be approximated by a percentage of impervious surface was greatest in Macao, followed by Dongguan, Shenzhen, Zhongshan, Foshan, Zhuhai, Guangzhou, Hong Kong, Jiangmen, Huizhou, and Zhaoqing. The last four cities had the percentage of impervious surface smaller than that in the GBA, suggesting that they had more natural land sources. The percentages of different land cover types in the GBA and different cities are presented in Figure 5 to reveal land cover composition. Results indicated that forest, cropland, and impervious surface were the three major land cover types in the GBA. As for the eleven cities, forests were the primary land cover in Guangzhou, Huizhou, Shenzhen, Zhuhai, Jiangmen, Zhaoqing, and Hong Kong; and the other four cities (i.e., Dongguan, Foshan, Zhongshan, and Macao) had the largest proportion in impervious surface. Land urbanization that can be approximated by a percentage of impervious surface was greatest in Macao, followed by Dongguan, Shenzhen, Zhongshan, Foshan, Zhuhai, Guangzhou, Hong Kong, Jiangmen, Huizhou, and Zhaoqing. The last four cities had the percentage of impervious surface smaller than that in the GBA, suggesting that they had more natural land sources.
Overall, the GBA was a megacity region with a high urbanization level and abundant natural landscape. In addition, the composition and distribution of the different land cover types were varied across the GBA and among different cities. Therefore, it is important to conduct an analysis focusing on the spatial pattern of the ecosystem services. Overall, the GBA was a megacity region with a high urbanization level and abund natural landscape. In addition, the composition and distribution of the different l cover types were varied across the GBA and among different cities. Therefore, it is portant to conduct an analysis focusing on the spatial pattern of the ecosystem servic

Ecosystem Services of the GBA and Different Cities
The capacity scores of ecosystem services supply, demand, and budget in the G are presented in Table 2. The total budget in the GBA was 184,467,815, indicating supplied ecosystem services were sufficient to meet demand capacity. The GBA had highest supply capacity in regulating services, followed by provisioning and cultural vices. There was the highest demand capacity in provisioning services for the GBA, lowed by regulating and cultural ones. Great differences in ecosystem services among different cities were observed shown in Figure 6. Although budgets were positive in most cities, there are several ex tions, for example, all the services categories in Macao, provisioning service Zhongshan, Dongguan, Shenzhen, and Foshan. The total budgets for Macao, Zhongs Dongguan, and Foshan were also negative. Among the 11 cities, the supply capacit regulating, provisioning, and cultural services was the largest in Zhaoqing, followed Huizhou, Jiangmen, and Dongguan. The top four cities that had higher demand capa were Guangzhou, Huizhou, Jiangmen, and Foshan. Focusing on the three services c gories, the supply capacity for the different cities tended to share a similar pattern: re lating > provisioning > cultural services. The demand capacity also had the same desce ing order as provisioning > regulating > cultural services for all the 11 cities.

Ecosystem Services of the GBA and Different Cities
The capacity scores of ecosystem services supply, demand, and budget in the GBA are presented in Table 2. The total budget in the GBA was 184,467,815, indicating that supplied ecosystem services were sufficient to meet demand capacity. The GBA had the highest supply capacity in regulating services, followed by provisioning and cultural services. There was the highest demand capacity in provisioning services for the GBA, followed by regulating and cultural ones. Great differences in ecosystem services among different cities were observed, as shown in Figure 6. Although budgets were positive in most cities, there are several exceptions, for example, all the services categories in Macao, provisioning services in Zhongshan, Dongguan, Shenzhen, and Foshan. The total budgets for Macao, Zhongshan, Dongguan, and Foshan were also negative. Among the 11 cities, the supply capacity of regulating, provisioning, and cultural services was the largest in Zhaoqing, followed by Huizhou, Jiangmen, and Dongguan. The top four cities that had higher demand capacity were Guangzhou, Huizhou, Jiangmen, and Foshan. Focusing on the three services categories, the supply capacity for the different cities tended to share a similar pattern: regulating > provisioning > cultural services. The demand capacity also had the same descending order as provisioning > regulating > cultural services for all the 11 cities.

Ecosystem Services along the Urban-Rural Gradient
The mean value of different ecosystem services capacity along the urban-rural gradient in the GBA is shown in Figure 7. Supply and demand capacity showed an increasing and decreasing trend from urban center to rural areas, respectively, which is attributed to the fact that the urban core in the GBA is agglomerated by the impervious surface of several cities and surrounded by natural landscape such as forests. The peak and valley points for supply capacity of the three services categories turn to be the 22nd and 3rd zones, Sustainability 2021, 13, 7044 9 of 16 which are the opposite for demand capacity. Therefore, the 22nd and 3rd zones have the largest and smallest total budgets, which correspond to regions with the highest percentage of the natural landscape and impervious surface, respectively.

Ecosystem Services along the Urban-Rural Gradient
The mean value of different ecosystem services capacity along the urban-rural dient in the GBA is shown in Figure 7. Supply and demand capacity showed an increa and decreasing trend from urban center to rural areas, respectively, which is attribut the fact that the urban core in the GBA is agglomerated by the impervious surface of eral cities and surrounded by natural landscape such as forests. The peak and v points for supply capacity of the three services categories turn to be the 22nd and zones, which are the opposite for demand capacity. Therefore, the 22nd and 3rd z have the largest and smallest total budgets, which correspond to regions with the hig percentage of the natural landscape and impervious surface, respectively.

Ecosystem Services along the Urban-Rural Gradient
The mean value of different ecosystem services capacity along the urban-rural gradient in the GBA is shown in Figure 7. Supply and demand capacity showed an increasing and decreasing trend from urban center to rural areas, respectively, which is attributed to the fact that the urban core in the GBA is agglomerated by the impervious surface of several cities and surrounded by natural landscape such as forests. The peak and valley points for supply capacity of the three services categories turn to be the 22nd and 3rd zones, which are the opposite for demand capacity. Therefore, the 22nd and 3rd zones have the largest and smallest total budgets, which correspond to regions with the highest percentage of the natural landscape and impervious surface, respectively. Furthermore, the mean value of different ecosystem services capacity along the urban-rural gradient for the 11 cities is presented in Figure 8. It should be noted that the construction of buffer zones for each city does exclude regions of other cities to assess Furthermore, the mean value of different ecosystem services capacity along the urbanrural gradient for the 11 cities is presented in Figure 8. It should be noted that the construction of buffer zones for each city does exclude regions of other cities to assess spatial neighboring effects. The first zone that covers other cities is marked with a vertical dotted line, and results excluding other cities are presented by dotted lines as a comparison group (Supply_com, Demand_com, and Budgets_com in Figure 8). The vast majority of buffer zones around urban centers covered other cities, for example, 18 out of 21 for Guangzhou and 19 out of 20 for Zhuhai, since the urban centers are close to each other with the mean nearest neighbor distance about 60 km. In view of this, relatively greater accessibility of ecosystem services for urban centers can be provided by urban areas in other cities rather than its own rural area. Therefore, it is vital to exploit the spatial neighboring effects. In contrast to comparison groups, neighboring cities generally had a negative effect on the ecosystem service according to the smaller budgets. However, for Shenzhen zones within 20-30 km from an urban center, a positive neighboring effect was observed. The general increasing patterns for budgets from urban centers to rural areas were influenced in some cities. There were two types of exceptions: fluctuation (i.e., Huizhou and Zhuhai) and decrease (i.e., Jiangmen, Zhaoqing, and Hong Kong), which may be attributed to the existence of multiple centers. In addition, the three cities (i.e., Jiangmen, Zhaoqing, and Hong Kong) provided supply that exceeds demand in all the zones, suggesting the requirement of ecological benefits for humans in both urban and rural areas are met. It is interesting to note that Zhongshan and Macao were balanced or undersupplied in all the zones. The other cities had an undersupply of ecosystem services in the urban core area, and other zones were characterized by ecosystem services' supplies exceeding their demands.
Sustainability 2021, 13, x FOR PEER REVIEW 10 of 16 spatial neighboring effects. The first zone that covers other cities is marked with a vertical dotted line, and results excluding other cities are presented by dotted lines as a comparison group (Supply_com, Demand_com, and Budgets_com in Figure 8). The vast majority of buffer zones around urban centers covered other cities, for example, 18 out of 21 for Guangzhou and 19 out of 20 for Zhuhai, since the urban centers are close to each other with the mean nearest neighbor distance about 60 km. In view of this, relatively greater accessibility of ecosystem services for urban centers can be provided by urban areas in other cities rather than its own rural area. Therefore, it is vital to exploit the spatial neighboring effects. In contrast to comparison groups, neighboring cities generally had a negative effect on the ecosystem service according to the smaller budgets. However, for Shenzhen zones within 20-30 km from an urban center, a positive neighboring effect was observed. The general increasing patterns for budgets from urban centers to rural areas were influenced in some cities. There were two types of exceptions: fluctuation (i.e., Huizhou and Zhuhai) and decrease (i.e., Jiangmen, Zhaoqing, and Hong Kong), which may be attributed to the existence of multiple centers. In addition, the three cities (i.e., Jiangmen, Zhaoqing, and Hong Kong) provided supply that exceeds demand in all the zones, suggesting the requirement of ecological benefits for humans in both urban and rural areas are met. It is interesting to note that Zhongshan and Macao were balanced or undersupplied in all the zones. The other cities had an undersupply of ecosystem services in the urban core area, and other zones were characterized by ecosystem services' supplies exceeding their demands.

Ecosystem Services Balance at Local Scale
Ecosystem services balance at multiple scales was analyzed by presenting the percentage of different groups of balance (i.e., negative values indicate demand exceeds supply/undersupply, and positive values indicate supply exceeds demand/oversupply) in Figure 9. The following groups were considered: <−4 = very high undersupply, (−4, −3) = high undersupply, (−3, −2) = medium undersupply, (−2, −1) = low undersupply, (−1, 0) = very low undersupply, (0, 1) = very low oversupply, (1, 2) = low oversupply, (2, 3) = medium oversupply, (3,4) = high oversupply, >4 = very high oversupply. The influence of the neighborhood size on the percentages of different balance groups is variant for different land cover types. Several general observations can be drawn from Figure 9. Over 90% of the impervious surface had very high undersupply, and over 90% of forests had a very high oversupply. The proportion of oversupply was over 80% on the scale of 1000 m and increased over 90% on smaller scales (i.e., 500 m and 200 m). Cropland and grassland showed a descending trend for percentages from very low to very high oversupply groups. Their undersupply also took up about 5% and 20%, respectively. It is interesting to note that the inclusion of neighborhood increased demand scores in a small proportion of impervious areas since the impervious surface had a very low budget score and its density was high in the urban area. Another land cover type (i.e., bare land) also had a negative balance score and showed oversupply in a very large percentage region when neighborhood landscape was included.

Ecosystem Services Balance at Local Scale
Ecosystem services balance at multiple scales was analyzed by presenting the percentage of different groups of balance (i.e., negative values indicate demand exceeds supply/undersupply, and positive values indicate supply exceeds demand/oversupply) in Figure 9. The following groups were considered: <−4 = very high undersupply, (−4, −3) = high undersupply, (−3, −2) = medium undersupply, (−2, −1) = low undersupply, (−1, 0) = very low undersupply, (0, 1) = very low oversupply, (1, 2) = low oversupply, (2, 3) = medium oversupply, (3,4) = high oversupply, >4 = very high oversupply. The influence of the neighborhood size on the percentages of different balance groups is variant for different land cover types. Several general observations can be drawn from Figure 9. Over 90% of the impervious surface had very high undersupply, and over 90% of forests had a very high oversupply. The proportion of oversupply was over 80% on the scale of 1000 m and increased over 90% on smaller scales (i.e., 500 m and 200 m). Cropland and grassland showed a descending trend for percentages from very low to very high oversupply groups. Their undersupply also took up about 5% and 20%, respectively. It is interesting to note that the inclusion of neighborhood increased demand scores in a small proportion of impervious areas since the impervious surface had a very low budget score and its density was high in the urban area. Another land cover type (i.e., bare land) also had a negative balance score and showed oversupply in a very large percentage region when neighborhood landscape was included.

Discussion
A series of policies have been launched in the GBA to promote ecological land protection (e.g., National Forest City Group Construction in the Pearl River Delta (2016-2025)). Benefiting from the impact of National Forest City Group Construction, the forest land reached 62.5%, according to our land cover map in 2020. As reported by a recent study, ecological land shifted from decline and fragmentation during 1990-2010 to growth

Discussion
A series of policies have been launched in the GBA to promote ecological land protection (e.g., National Forest City Group Construction in the Pearl River Delta (2016-2025)). Benefiting from the impact of National Forest City Group Construction, the forest land reached 62.5%, according to our land cover map in 2020. As reported by a recent study, ecological land shifted from decline and fragmentation during 1990-2010 to growth and integration during 2010-2019, and 659 km 2 of non-ecological land was converted to ecological land, which may be in response to the "ecological red line" policy [33]. Since forest provides the largest budget capacity among the seven land cover types, the GBA can obtain an oversupply of ecosystem services with rich forest resources. Although there were very large budgets for regulating, provisioning, and cultural services in the GBA, city-scale results indicated a quietly unbalanced supply-demand mismatch. Specifically, only Zhaoqing, Huizhou, and Jiangmen had relatively greater supply capacity than demand capacity, and the supply capacity of provisioning services in Zhongshan, Dongguan, Shenzhen, and Foshan could not meet the demand. Meanwhile, Macao faces serious ecological security as all three general services categories were negative. Since cropland, forest, shrubland, and water can supply more provisioning services, it is very necessary to protect or even increase these areas in the above-mentioned four cities. Corresponding city-scale policies have been considered. Shenzhen City issued a "Management stipulation of the basic ecological line" that forbids any construction projects within the basic ecological line [34]. Land policy in Hong Kong emphasizes ecological protection, and its Town Planning Ordinance constitutes the legal basis for the conservation or protection of country parks, green belts [35], making it a highly urbanized city with favorable ecological conditions. How to further improve the effectiveness of ecological land protection policies has become a major challenge.
According to the Outline Development Plan for the Guangdong-Hong Kong-Macao Greater Bay Area, the development goal for Huizhou is to turn into a green eco-tourism city and to take ecological responsibility for the Greater Bay Area [36]. The ecosystem service supply capacity of Huizhou ranked second among the 11 cities, which followed behind Zhaoqing. However, Huizhou has a unique location (i.e., adjacent to Guangzhou, Dongguan, and Shenzhen (they are developed cities with dense population density and intensive industry)). As indicated by the urban-rural gradient analysis, interregional cooperation of Shenzhen and other neighboring cities lead to improved budgets in the urban zones. Therefore, interregional cooperation of the multiple cities in the GBA from the perspective of ecological conservation is encouraged. With the growing coordination development of various aspects like economic development, public services, infrastructure construction, and environmental protection in the GBA [37,38], spatial fairness of ecosystem services at the city and local scale should be improved by putting forward multi-city ecological protection and coordinated development policies. Our spatial pattern analysis can provide essential information for further regional development plans without administrative constraints.
The spatial pattern analysis in this study is only concerned about the ecosystem in each administrative district scale (i.e., the GBA and 11 cities), along the urban-rural gradient, and at a local scale. To deal with it, most manual spatial pattern analyses (e.g., zone buffers and local neighborhood analysis) were conducted. Recently, some more advanced spatial pattern mining techniques (e.g., spatial hotspot detection, colocation detection, and outlier detection [39], and even spatio-temporal pattern analysis [40]) have been utilized in other studies and show promising results. Therefore, more in-depth spatial and spatio-temporal pattern analysis work can be considered in further research.

Conclusions
With the support of Sentinel-1 SAR and Sentinel-2 optical remote sensing data, we proposed a deep learning-based land cover mapping framework to generate a 10 m spatial resolution map. Based on the land cover map of the GBA in 2020, the fine-scale spatial pattern of ecosystem services was analyzed. Specifically, total ecosystem services supply, demand, and budgets in the GBA and the 11 cities are presented. Ecosystem services along the urban-rural gradient were analyzed to reveal their trend from urban centers to rural areas. Finally, ecosystem services balance assessments at the local scale were performed to identify the pixel-level status of supply and demand balance. The major study results showed as follows: (1) Forest, cropland, and impervious surface were the three major land cover types in the GBA. Forest was the primary land cover in Guangzhou, Huizhou, Shenzhen, Zhuhai, Jiangmen, Zhaoqing, and Hong Kong, and the impervious surface was the main land cover in the other four cities. (2) Although ecosystem services in the GBA were sufficient to meet their demand, there was undersupply for all the three general services in Macao, provision services in Zhongshan, Dongguan, Shenzhen, and Foshan. (3) Along the urban-rural gradient in the GBA, supply and demand capacity showed an increasing and decreasing trend since the urban core in the GBA is agglomerated by the impervious surface of several cities and is surrounded by natural landscape such as forests. As for the city-level urban-rural gradient analysis, except for the general increasing pattern for budgets from the urban centers to rural areas, Huizhou and Zhuhai showed a fluctuation pattern, while Jiangmen, Zhaoqing, and Hong Kong presented a decreasing pattern, which may be attributed to the existence of multiple urban centers. (4) The inclusion of neighborhood landscape increased demand scores in a small proportion of impervious areas, and bare land with a negative balance score showed an oversupply in a very large percentage region when neighborhood landscape was included.
To realize the sustainable development of the GBA, on the one hand, local governments should establish policies that strictly protect ecological land (i.e., forests, shrubland, water, and cropland); on the other hand, interregional cooperation of the multiple cities in the GBA from the perspective of ecological conservation is encouraged.