Spatial Correlation between Ecosystem Services and Human Disturbances: A Case Study of the Guangdong–Hong Kong–Macao Greater Bay Area, China

Abstract

Exploring the spatial relationship between ecosystem services (ES) and human disturbance intensity (HDI) is vital for maintaining regional ecological security. This study aims to explore the spatial correlation between ES and HDI in the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) and provide meaningful implications for coastal ecological planning. Multi-source remote sensing data, remote sensing software, and geographic information system provided initial data and technical support for this research. We integrated four human pressures (population, land-use, traffic, and energy) to map the HDI in the GBA for 2018. Coastal ES were comprehensively considered and spatially visualized by extracting the ES sources. The geographically weighted Pearson correlation coefficient and bivariate local Moran were used to quantitatively reflect and spatially visualize the detailed relationship between ES and HDI. Our study presents several key findings. First, the central and southern parts of the GBA are under strong HDI, dominated by a dense population and intense land utilization. Second, the kernel density of ES sources can better manifest the spatial distribution of ES objectively in comparison to the traditional model calculation. Provisioning services mainly originate from the periphery of the central cities; cultural services are highly concentrated in the heartland of the GBA; and regulating and maintenance services have high density in the outermost regions. Third, ES and HDI have a significant correlation, and the geographically weighted Pearson correlation coefficient and local indicator of spatial association cluster maps illustrate that unlike the global findings, the local correlation is spatially nonstationary as the local scale is affected by specific human activities, natural conditions, regional development, and other local factors. Four, high-capacity regions of ES provision are mainly under high HDI. Areas with high provisioning service values are mainly affected by population and traffic pressure, whereas regulating and maintenance services and cultural services are mainly dominated by high-density populations. Regulating and maintenance services are also affected by land-use pressure. We determine that human disturbance has negative spillover effects on ES, which should be the focus in regional ecological planning.

1. Introduction

Ecosystem services (ES) are defined as the benefits people directly or indirectly obtain from ecosystems; they are classified into three categories by the Common International Classification of Ecosystem Services (CICES): provisioning, regulating, and maintenance, and cultural services [1]. These services are a necessary basis for socioeconomic development; ES abundance is also a resilience descriptor for ecological security [2]. Intensified human disturbance and the tight-coupled human–environment system have dramatic impacts on the quality and quantity of ES [3,4,5,6,7]. Coastal areas are economically important and provide abundant ES. However, coastal regions are ecologically fragile and sensitive to environmental change [8]. The vulnerability of coastal ecosystems has been increasing as a result of human perturbations that compromise the ability to provide multiple ecosystem services. Human pressures caused by natural resource exploitation, infrastructure construction, and sea reclamation can exert irreversible effects on coastal ecosystems [9,10,11]. Therefore, exploring the spatial relationship between ES and human disturbance intensity (HDI) is fundamental to making informed decisions regarding regional ecological security and alleviating the contradiction between ES supply and human demand.

The dynamic mechanism between human activities and ES has been a topic of considerable interest in studies on human–nature interaction systems. Many studies have been conducted to construct a globally standardized measure of the human footprint on natural ecosystems. Sanderson et al. (2002) mapped the intensity of human influence over the entire land surface based on remote sensing data from the early 1990s [12]. It is the most complete and with the highest resolution global dataset of accumulated human pressure on the environment. Halpern (2008) developed an ecosystem-specific and multiscale spatial model to synthesize 17 global datasets of anthropogenic drivers of ecological change for 20 marine ecosystems [13]. Venter et al. (2016) used 1 km² resolution remote sensing data from 1993 to 2009 of infrastructure, land cover and human access to natural areas to construct a globally standardized measure of the cumulative human footprint on the terrestrial environment [14]. Tapia-Armijos et al. (2017) tailored the Human Footprint Index at the global level to evaluate the spatial and temporal patterns of human pressure in the protected areas of South Ecuador [15]. These studies indicate that human demands on natural systems are accelerating and could be undermining the stability of ecosystems, and that most regions are strongly affected by multiple artificial drivers. They also suggest that remote sensing can play a significant role in human–nature interaction research for important initial data and technical support.

The main human disturbance factors vary between spatial scales, and differ by geographical location; furthermore, global datasets and measurements often overlook information that may be unique to specific localities. Rogers et al. (2010) incorporated human-related threats obtained by remote sensing data processing with ES and biological measures to identify priorities within a developing protected area network [16]. Li (2018) used human footprint scores processed by Landsat datasets to reveal the spatial relationships between the human footprint and the environmental quality of nature reserves in Tibet [17]. Mustafa and Blaine (2020) examined the variation in habitat quality associated with HDI [18]. Although these studies attempted to explore the interaction between ES and HDI in local areas, there are still some factors that remain unexplored (like the clustering patterns between ES and HDI), especially at the local scale. Peng et al. (2017) used linear regression to estimate the linear relationship between the total ES value and three indicators of urbanization [19]. Nevertheless, linear regression and the global Moran method for correlation analysis can merely analyze the correlation quantitatively, and the correlation cannot be visualized in space [20,21].

The aforementioned studies were conducted in inland nature reserves and habitats; research on sensitive coastal areas with important ecological functions is relatively scarce. Compared with terrestrial ecosystems, coastal areas have more complex environmental conditions and ecosystem components. Existing studies have mainly estimated and mapped the value of ES [7,22,23], which seriously underestimate the actual value of ES in coastal areas. Therefore, determining ways to effectively and comprehensively identify the types and distribution characteristics of ES in coastal urban areas is of great significance for ecological planning, and ecosystem function restoration and improvement in coastal zones. Hu et al. (2019) used the hot/cool-spot mapping method to assess the patterns of ES values in the Pearl River Delta, and discussed the response mechanism of ES value to land-use/cover changes [24]. However, there have been no studies on human-nature interactions in the coastal region of south China.

This study aims to explore the spatial relationship between ES and HDI in the Guangdong–Hong Kong–Macao Greater Bay Area (GBA) in 2018 and provide foundations for mitigating environmental damage in sensitive or ecologically valuable coastal areas. The objectives of this study are to: (1) map the cumulative HDI by incorporating four types of human pressures based on multi-source remote sensing data; (2) depict the spatial pattern of ES according to the coastal ES classification system; (3) analyze the spatial correlations and clustering patterns of ES and HDI; and (4) identify the implications for coastal ecological planning.

The GBA is the fourth largest bay area in the world, and it is one of the most economically vibrant regions in China. In recent decades, large-scale land exploitation and reclamation have been carried out, which have resulted in irreversible adverse effects on the coastal ecological environment [25]. To ensure the harmonious and sustainable development of human society and environmental preservation in the GBA, it is important to explore the spatial correlation between HDI and ES. In this study, we mapped the cumulative HDI by incorporating four main types of human pressures in the GBA (population, land-use, traffic, and energy). Pressures caused by navigable waterways were included to ensure the HDI calculation reflected the specific context of coastal zones. The spatial pattern of ES was depicted according to the coastal ES classification system; the multiple types of ES provided by coastal ecosystems were fully considered. The density analysis of points of interest could reflect urban economic and population hotspots [26]. Therefore, we extracted the ES source points to visualize the spatial distribution of ES hotspots and calculated the kernel density for the division of regional ES capacity. The kernel density of ES sources surpasses the subjectivity of traditional modeling methods and manifests the spatial distribution of ES more objectively [27]. The spatial correlations and clustering patterns of ES and HDI were analyzed by applying the geographically weighted Pearson correlation coefficient and the local spatial autocorrelation, which can spatially visualize and quantitatively reflect the correlation between the two variables in the local area.

2. Materials and Methods

2.1. Study Area

The GBA is located in southern China (21°25′ N–24°30′ N, 111°12′ E–115°35′ E) and covers a total area of 56,000 km² (Figure 1). It consists of the Pearl River Delta and two special administrative regions (Hong Kong, and Macao); the region plays a significant strategic role in the land-based “Silk Road Economic Belt,” and “the 21st Century Maritime Silk Road” [28]. The GBA is characterized by subtropical and tropical monsoon climates; it records an annual average precipitation ranging from 1300 mm to 2500 mm, an average temperature of 22.3 °C, and exhibits two primary seasons: a high-temperature summer with more rain, and a milder winter, with less rain. The altitude is higher in the north, and lower in the south. The mountains are mainly distributed in the northeast and northwest regions; the central and coastal zones are dominated by plains [29].

2.2. Overview of the Methodological Steps

In this study, four main types of human disturbance were identified and superposed to obtain the cumulative HDI. Three types of ES were identified based on the coastal ES classification system and were extracted via the points of interest extractor. The spatial correlations and clustering patterns of ES and HDI were analyzed using the geographically weighted Pearson correlation coefficient in R Programming language and bivariate local Moran in GeoDa software. The flow chart of the entire process is shown in Figure 2.

2.3. Human Disturbance Intensity

The human population of the GBA far exceeds that of the other three world-class bay areas (New York Bay Area, San Francisco Bay Area, and Tokyo Bay Area); in 2018, the GBA recorded a population density of 881.17 people/km² [29,30]. Population growth has increased the burden on the local ecological environment.

A large number of reclamation projects have been conducted in coastal areas and the offshore islands of the GBA. The built-up land covers an area of 6769.17 km², which accounts for 12.18% of the total land area of the GBA. High-intensity human activities have led to changes in land cover to varying degrees and have greatly affected the ecosystem [31,32]. We preprocessed the remote sensing images (Landsat 8 OLI/TIRS) using ENVI 5.3, including performing geometric, radiometric, atmospheric corrections, image cropping, and mosaics, and then classified the preprocessed images through support vector machine supervised classification to obtain the land-use data. We combined Google Earth data with land-use data to establish a confusion matrix to verify the accuracy of the classification results. The verification results show a Kappa coefficient >0.85, which means the interpretation accuracy meets the use requirements of medium-resolution remote sensing images.

Dense transportation networks negatively affect the integrity of the ecological landscape, aggravate landscape fragmentation, and weaken landscape connectivity [33]. In 2018, the GBA had 4300 km of highways, and the road network density was much higher than that of New York, Tokyo, London, the Yangtze River Delta, the Beijing–Tianjin–Hebei region, and other domestic and foreign city clusters [28,29,30]. In addition, the degree of connectivity of the night-time lights in the GBA is significantly higher than in other urban agglomeration.

Satellite images capture visible lights at night and indicate the degree of development of an area [34]. There is a strong correlation between power consumption intensity and night-time light disturbance [35,36].

Considering the aforementioned observation, we took the four types of human disturbance as leading contributors to cumulative HDI; the data sources are shown in

Table 1. Multi-source remote sensing data used in the study.

Data Name	Data Type	Spatial Resolution/Scale	Source
Population density data	Remote sensing data	1 km grid	LandScan Global Population Database, Department of Energy, Oak Ridge National Laboratory [37]
Landsat 8 OLI/TIRS images	Remote sensing data	30 m	United States Geological Survey [38]
Road network vector data	Vector data	1:1,000,000	China National Catalogue Service for Geographic Information [39]
VIIRS night-time lights image	Remote sensing data	500 m	Google Earth Engine (GEE, https://code.earthengine.google.com accessed on 18 March 2021)

. We used a 3 km × 3 km vector grid as the spatial unit for analysis to extract the mean value of multi-source data in ArcGIS. Then, the disturbance intensity scores were calculated according to the corresponding formula.

2.3.1. Disturbance Intensity of Population

We assigned a disturbance intensity score of 10 for grids with more than 1000 people/km²; this was based on the assumption that the influence of the human population approached an asymptote at 1000 people/km² [3]. Owing to the higher population density in developed coastal areas, the asymptote was set to 1000 rather than 10 inhabitants/km² as was used in global datasets. For pixels with population density below 1000 people/km², the population disturbance score $\left(PD{I}_S\right)$ was determined using the calculation method as in global datasets [14].

$$PD{I}_s=\left\{\begin{array}{c} 10, Pd\ge 1000 \text{people} / {\mathrm{km}}^2\\ 3.333\times \log \left(Pd+1\right), Pd<1000 \text{people} / {\mathrm{km}}^2\end{array}\right.$$

(1)

The value of $PD{I}_S$ is in the range (0, 10]; the higher the score, the stronger is the disturbance intensity.

$$\mathrm{Nor}PD{I}_S\left(i\right)=\frac{PD{I}_S\left(i\right)-\min PD{I}_s}{\max PD{I}_s-\min PD{I}_s},$$

(2)

where Pd is the population density, and where $\mathrm{Nor}PD{I}_S\left(i\right)$ is the normalized $PD{I}_S$ of each grid, and ranges from 0 to 1. $PD{I}_S\left(i\right)$ is the $PD{I}_S$ of grid i, $\max PD{I}_s$ is the maximum of $PD{I}_s$ in all grids, and $\min PD{I}_s$ is the minimum of $PD{I}_s$ in all grids.

2.3.2. Disturbance Intensity of Land-Use

The areas of each land-use type in each grid were calculated using the research method developed by Brentrup et al. (2002) [40]. The calculation formula for the land-use intensity is as follows:

$$LD{I}_s={{\displaystyle \sum }}_{i=1}^n(P_i\times S_i),$$

(3)

where $LD{I}_S$ is the land-use intensity score, n is the total number of land-use types included in a grid, P_i is the percentage of land-use type i to the total area of a grid, and S_i is the pressure score of the land-use type i [41] (

Table 2. Land-use types and land-use intensity scores assigned to each type.

Land-Use Type	Description	Disturbance Intensity	Score
Built-up land	Land for urban uses, including factories, quarries, mining, transportation facilities, and airport	excessive	10
Rural settlement	Land used for settlements in villages	very strong	8
Cultivated land	Land mainly for growing crops	strong	6
Garden land	Unformed forest afforestation land, nursery, and various garden land		6
Inland water	Natural inland waters, and land for water conservancy facilities, including rivers, lakes, reservoirs, ponds, etc.	moderate	4
Grassland	Area dominated by herbaceous plants of natural growth and artificial planting		4
Woodland	Tracts of natural forest, secondary forest, and artificial forest	Weak	2
Coastal wetland	Depressions in the coastal area that has been in a state of standing water or semi-water for a long time, including beaches, swamps, mangroves, etc.	Slight	1
Unused land	Not put into practical use or is difficult to use, including sandy land, saline land, swampland, bare soil, bare rock, and others	Almost no disturbance	0

). We assigned the maximum score of 10 to built-up land, which is the land-use type most influenced by human activities. Lower scores were allocated to woodland, coastal wetland, and unused land. We normalized the land-use intensity score in each grid to obtain $\mathrm{Nor}LD{I}_s\left(i\right)$.

2.3.3. Disturbance Intensity of Transportation

Buffer zones of different radii were created in ArcGIS to depict different levels of traffic pressure in the GBA. We assigned influence scores from 0–10 to each transport line falling within different buffer ranges (

Table 3. Traffic disturbance intensity scores.

Type	0–1 km	1–5 km	5–10 km	10–15 km
Expressway	10	8	7	5
First-grade highway	10	8	4	2
Secondary highway	8	6	2	0
Tertiary highway	6	4	1	0
Fourth-class highway	4	2	0	0
Substandard way	2	1	0	0
Rural road	1	0	0	0
Railway	9	7	5	3
Navigable waterway	7	5	1	0

) [14,18]. The higher the score, the stronger is the disturbance intensity. Then, the total traffic disturbance intensity score was obtained by adding the scores of each transport line in each grid. Similarly, we standardized the traffic pressure in each grid to obtain $\mathrm{Nor}TD{I}_S\left(i\right)$.

Navigable waterways are also significant transportation corridors in the GBA; in 2018, the total mileage of navigable waterways in Guangdong Province reached 16,853 km. The cargo throughput of GBA port group in 2018 was 1.73 billion tons [29]. Therefore, we considered the disturbance of the navigable waterways in our study and assigned pressure scores to the buffer zones [14]. We only considered the main navigable waterways in the Pearl River Delta region because the waterway transportation is almost entirely along the coast in Hong Kong and Macao.

2.3.4. Disturbance Intensity of Energy Consumption

The data of night-time lights were converted into the energy consumption intensity. The inter-calibrated digital number of night-time lights was rescaled using an equal quantile approach in ArcGIS on a 0–10 scale [14]. The higher the score, the stronger is the disturbance intensity. We normalized the energy consumption intensity in each grid to obtain $\mathrm{Nor}ED{I}_S\left(i\right)$.

2.3.5. Cumulative Human Disturbance Intensity

Our research on the various human disturbance factors of the GBA did not involve specific weighting; therefore, to maintain research neutrality, the default weights of the four human disturbance indices were kept the same. The cumulative HDI score of grid i is the summation of the four human pressure scores. The expression is as follows:

$$HDI\left(i\right)=\mathrm{Nor}PD{I}_S\left(i\right)+\mathrm{Nor}LD{I}_S\left(i\right)+\mathrm{Nor}TD{I}_S\left(i\right)+\mathrm{Nor}ED{I}_S\left(i\right).$$

(4)

The cumulative HDI was classified using Jenks natural breaks method into five HDI classes: very low, low, moderate, high, and very high.

2.4. Coastal Ecosystem Services System

Coastal ecosystems are composed of terrestrial, freshwater, wetland, and nearshore ecosystems. The coastal ES classification system was constructed based on the CICES [1]. The three categories and 12 subclasses of ES were identified from the aforementioned four types of ecosystems, and the ES sources were selected with reference to Li and Wang (2019) [42] (

Table 4. Ecosystem service classification system and the ecosystem service sources of the GBA.

Main Category	Subclass	Specific Content	ES Source
Provisioning services	Nutrition and essential substances	Terrestrial and aquatic animal and plant food, drinking water	Orchard, breeding base, reservoir
	Raw material	Biological and non-biological materials (e.g., wood, rubber)	Forestry center
	Energy sources	Renewable biofuels and non-bioenergy (such as minerals, geothermal resources, hydropower, solar energy)	Mine (including geothermal energy), power station
	Other	Carrier of other service functions (such as transport carrier)	Port wharf
Regulating and maintenance services	Hydrology regulation and water purification	Intercept, absorb, and store precipitation, regulate runoff; filter and decompose impurities and harmful chemicals, and promote the self-purification of the water	Natural habitat
	Climate regulation and air quality maintenance	Promote the self-purification of the air, adjust the climate, and provide a climate suitable for human survival	Natural habitat
	Soil conservation, wind break, and sand fixation	Retain soil and slow down erosion; intercept and decompose organic matter, provide fertile land resources, and soil self-repair ability	Natural habitat
	Habitat maintenance and biodiversity conservation	Provide a stable habitat environment for organisms, carry out biological gene preservation, population maintenance, and new species breeding, and protect biodiversity	Nature reserve
	Flood and tidewater control and shoreline stabilization	Confers resistance to disasters, adapt to disturbances, and maintain ecosystem stability	Mangrove, coral reef
Cultural services	Landscape	Aesthetic appreciation, recreation, and leisure travel	Observation deck, bathing beach, forest park, zoo
	Cultural carrier	Spirit and culture contained in the ecosystem, such as history and religion
	Scientific research and education	The object of scientific research and education provides opportunities to understand, observe, and explore the ecosystem

).

Provisioning services refer to the capacity of the environment to provide for the basic material needs of human survival and development, such as food, raw materials, and energy. Regulating and maintenance services refer to the role of the ecosystem in regulating and restoring the ecological security of the atmospheric, water, and soil environments through a series of ecological processes. Coastal wetlands (mangroves and coral reefs) can enhance the resilience of coastal zones under conditions of global climate change. Nature reserves are of great significance to protecting wildlife, forest ecology, inland wetlands, and ocean coasts [43]. Cultural services refer to the function of the ecological system as a carrier of esthetic value and historical culture, a place for scientific research, and education. Natural attractions, forest parks, and zoos have high natural ornamental value and scientific research and educational value [44]. Coastal bathing beaches provide residents with a place for leisure.

Points of interest cover the location and attribute information of varied geographic entities [45]. The ES sources were extracted from Amap using OSpider 2.1 on an administrative unit. The accuracy was calculated to be >85% by observing the consistency between the service functions of the sample points of ES sources and reference points from Open Street Map (https://www.openstreetmap.org/ accessed on 18 March 2021); therefore, the mapping results can accurately reflect the spatial distribution characteristics of ES in the GBA. The kernel density estimation method does not make any assumptions regarding the data distribution. It is a method used to study the characteristics of data distribution from the data sample itself [46]. It can be used to intuitively analyze the degree of aggregation of the ES spatial number distribution. We used the natural breaks (Jenks) classification method for the partitioning of kernel density to determine the capacity level of the ES in the GBA.

2.5. Spatial Correlation

Ecological and societal processes and patterns have uneven and complex spatial distributions due to differences in geographic environments. A geographically weighted model can be used to explore this characteristic among different variables [47,48,49]. We calculated the geographically weighted Pearson correlation coefficient [50] to explore the local correlation between ES and HDI. The formula is as follows:

$$GWPCC=\frac{{{\displaystyle \sum }}_{j=1}^n{w}_{ij}\left(x_j-{\overline{x}}_i\right)\left(y_j-{\overline{y}}_i\right)}{\sqrt{{{\displaystyle \sum }}_{j=1}^n{w}_{ij}{\left(x_j-{\overline{x}}_i\right)}^2}\sqrt{{{\displaystyle \sum }}_{j=1}^n{w}_{ij}{\left(y_j-{\overline{y}}_i\right)}^2}},$$

(5)

where GWPCC is the geographically weighted Pearson correlation coefficient of ES and HDI; $x_i$ and $y_i$ are the ES and HDI of the regression analysis point i, respectively; and $x_j$ and $y_j$ are the ES and HDI, respectively, of the neighboring point j. The number of nearest neighbors around observation point i is represented by n, and ${\overline{x}}_i$ and ${\overline{y}}_i$ are the geographically weighted mean value of ES and HDI of point i, respectively, which are calculated by the formulae:

$${\overline{x}}_i={{\displaystyle \sum }}_{j=1}^n{x}_j{w}_{ij}/{{\displaystyle \sum }}_{j=1}^n{w}_{ij},$$

(6)

$${\overline{y}}_i={{\displaystyle \sum }}_{j=1}^n{y}_j{w}_{ij}/{{\displaystyle \sum }}_{j=1}^n{w}_{ij},$$

(7)

$$w_{ij}=\left\{\begin{array}{c}{\left[1-{\left(d_{ij}/b\right)}^2\right]}^2, d_{ij}\le b\\ 0 , otherwise\end{array}\right.$$

(8)

where w_ij is the weight at point j, d_ij is the distance between points i and j, and b is the bandwidth of the kernel function. A smaller bandwidth improves the accuracy of the spatial change of the result, whereas a larger bandwidth makes the result appear closer to the global model [50]. By comparing the experimental results, it was found that setting the bandwidth to 0.03 produced results that are more in line with the spatial pattern of ES and HDI in the GBA.

The geographically weighted Pearson correlation coefficient ranged from −1 to 1. A positive value indicates a positive spatial correlation between ES and HDI, whereas a negative value indicates a negative spatial correlation. The greater the absolute value of the correlation coefficient, the stronger the spatial correlation between the ES and HDI.

2.6. Local Spatial Autocorrelation

Spatial autocorrelation analysis considers the relationship among spatially adjacent quantities. Bivariate local Moran analysis was performed in GeoDa V.1.12 to generate a local indicator of spatial association (LISA) cluster map, explore the spatial clustering patterns of ES and HDI, and verify their spatial autocorrelation. The clustering patterns were categorized into four types: high ES surrounded by high HDI (HH), high ES surrounded by low HDI (HL), low ES surrounded by high HDI (LH), and low ES surrounded by low HDI (LL).

3. Results

3.1. Spatial Patterns of Human Disturbance Intensity

Figure 3 demonstrates that human disturbances are unevenly distributed in the GBA. The strongest population disturbance is measured in the central cities (Hong Kong, Macao, Shenzhen, and Guangzhou); intermediate pressure is measured in Foshan, Zhongshan, Dongguan, and Zhuhai; and low pressure is measured in Jiangmen, Huizhou, and Zhaoqing (Figure 3a). The central cities are dominated by highly connected areas of built-up land, whereas cultivated land, garden land, and woodland are mainly distributed around the central cities and the periphery of the GBA. Built-up land exerts greater pressure on the ecological environment than other land types (Figure 3b). The traffic network system is radially distributed about the central cities. The intensity of traffic pressure in Guangzhou, Shenzhen, Foshan, and Huizhou is relatively high as compared to other cities (Figure 3c). The energy consumption of the central cities is significantly higher than that of the peripheral areas (Figure 3d). Central cities are international metropolises characterized by intense human activities and require huge amounts of energy to maintain their normal operations.

The center of the GBA is mainly dominated by moderate, high, and very high levels of HDI (Figure 4); these account for 22.33%, 15.92%, and 14.84% of the total intensity, respectively. The proportions of very low and low-intensity human disturbance are 21.95%, and 24.96%, respectively. The central and southern parts of the GBA are characterized by very high levels of human disturbance, especially in Guangzhou, Shenzhen, Zhuhai, Foshan, Dongguan, Zhongshan, Hongkong, and Macao, with mean HDI of 1.77, 2.38, 1.80, 2.00, 2.68, 2.35, 2.11, and 1.99, respectively. Some areas of Huizhou are characterized by medium and strong disturbance levels, whereas others exhibit only slight disturbance. Jiangmen and Zhaoqing have moderate and lower levels of disturbance, with mean HDI of 1.48, and 1.30, respectively. A very high-intensity disturbance in Guangzhou is distributed over its southern region. Most of the human disturbance in Hong Kong is of moderate or high intensity, but the average disturbance intensity score in its coastal segment is low. Levels of HDI are unevenly distributed over the GBA, as well as within the smaller administrative scale.

The contribution of each human pressure to the cumulative HDI was used to identify the leading contributors that affect the HDI level in the administrative area. Most cities are mainly affected by population density and intense land development (Figure 5). Compared with other disturbance factors, the population pressure of all municipal administrative regions is more prominent, and the population pressure of Guangzhou and Shenzhen is the most significant. The contribution of land-use pressure to the HDI of Guangzhou, Shenzhen, Zhuhai, Hong Kong, and Macao is relatively large at 30.83%, 27.42%, 33.03%, 30.24%, and 46.35%, respectively. The pressure of the traffic system provides greater contribution to the HDI in Guangzhou, Shenzhen, Foshan, Zhaoqing, Huizhou, Dongguan, and Zhongshan. Compared with other cities, Macao has higher energy consumption pressure.

3.2. Spatial Distribution Characteristics of Ecosystem Services

In the GBA, the delivery of food and nutrients is concentrated in the central built-up area and nearby vicinity with high population density (Guangzhou, Shenzhen, Zhuhai, Dongguan, and Foshan). The commercial forests and mineral resources are mostly distributed in Zhaoqing, Jiangmen, and Huizhou. Jiangmen and Dongguan provide the water and power resources necessary for human life and production (Figure 6a). The high- and medium-high-capacity areas of provisioning services account for 2.50% and 17.58% of the total area, respectively (Table 4). Medium capacity areas account for 26.30% of the total area, and low-medium and low-capacity areas account for 26.74%, and 26.88%, respectively. From a spatial perspective, the provisioning services originate from the periphery of the central cities (Figure 6b). Guangzhou, Shenzhen, Dongguan, Foshan, Jiangmen, and Huizhou have relatively high provisioning service value.

Mangrove ecosystems are mainly distributed in Guangzhou, Shenzhen, Zhuhai, Zhongshan, Jiangmen, and Huizhou (Figure 6c). The high- and medium-high- capacity areas of regulating and maintenance services account for 2.43% and 11.62% of the total area, respectively; the areas of medium capacity levels of regulating and maintenance services account for 19.36%. The low-medium and low-capacity areas account for 31.21%, and 35.38%, respectively. From a spatial perspective, the regulating and maintenance services show high density in the outermost periphery of the GBA (Figure 6d). The high and medium-high-capacity areas for regulating and maintenance services are distributed in the coastal areas of Guangzhou, Shenzhen, Zhuhai, Zhongshan, Zhaoqing, and Huizhou.

Leisure and sightseeing services are mostly distributed over the periphery of the central urban built-up areas and along the coastline (Zhuhai, Shenzhen, and Hong Kong). Cultural, scientific research, and education services are concentrated in the central cities that have relatively complete public service infrastructures (Figure 6e). The areas of high and medium-high-capacity levels of cultural services account for 2.29%, and 8.74% of the total area, respectively (

Table 5. Summary statistics of the capacity level of ecosystem services (ES).

Classification	Provisioning Services		Regulating and Maintenance Services		Cultural Services
	Area/km2	Proportion/%	Area/km2	Proportion/%	Area/km2	Proportion/%
low capacity	14,931	26.88%	19,647	35.38%	29,965	53.95
low-medium capacity	14,849	26.74%	17,334	31.21%	13,908	25.04
medium capacity	14,606	26.30%	10,753	19.36%	5542	9.98
medium-high capacity	9763	17.58%	6453	11.62%	4853	8.74
high capacity	1389	2.50%	1351	2.43%	1270	2.29

); cultural services are highly concentrated in the heartland of the GBA (Figure 6f). Medium capacity areas for cultural services account for 9.98% of the total area; low-medium and low-capacity cultural service areas account for 25.04%, and 53.95%, respectively.

3.3. Spatial Correlations between Ecosystem Services and Human Disturbance Intensity

The global bivariate Moran shows the strong negative correlation between provisioning services and HDI (bivariate Moran’s I: −0.2399, P < 0.01). The positive correlation is strongest between cultural services and HDI (bivariate Moran’s I: 0.4784, P < 0.01), followed by that between regulating and maintenance services and HDI (bivariate Moran’s I: 0.1665, P < 0.01).

The spatial patterns of the local correlations between ES and HDI are shown in Figure 7. There is significant negative correlation between provisioning services and HDI in Guangzhou’s central urban area, Shenzhen’s coastal zone, Foshan, Zhuhai, Hong Kong, Macao, and the north of Zhaoqing (Figure 7a,b). There is significant positive correlation in most of Zhongshan, Dongguan, Zhaoqing, Jiangmen, and Huizhou. There is significant negative correlation between RM and HDI in the coastal areas of Guangzhou, Shenzhen, Zhuhai, Zhongshan, and most areas of Zhaoqing, Jiangmen, and Huizhou (Figure 7c,d). Regulating and maintenance services degrade as human disturbance increase. There is significant positive correlation between CS and HDI in Guangzhou, Shenzhen, Zhuhai, Foshan, Huizhou, Hong Kong, and Macao (Figure 7e,f). The results show that highly developed, densely populated areas have more demand for cultural, educational, leisure, and entertainment services. Therefore, cultural services are highly exploited in these cities, and the density of cultural services increase as human disturbance increase.

A quantitative analysis was conducted on the area of space overlap for each capacity level of ES under different levels of HDI using the “tabulate area grid” statistical tool in ArcGIS 10.4 (Figure 8). The medium-high and high-capacity regions of ES in the GBA are mainly under high and very high HDI, respectively. A proportion of 46.56% medium capacity areas, 64.00% medium-high-capacity areas, and 95.25% high-capacity areas of provisioning services are affected by high and very high HDI. Over half of the medium capacity areas, medium-high, and high-capacity areas of regulating and maintenance services are under high and very high HDI. Among these, 62.36% are high-capacity areas affected by high to very high HDI. There are 45.42% low-medium regions, 70.57% medium capacity regions, 77.92% medium-high-capacity regions, and 75.02% high-capacity regions of cultural services, respectively, and these regions are influenced by high and very high HDI. Most of the low-capacity regions are under very low, low, or medium HDI.

3.4. Spatial Cluster Pattern of Ecosystem Services (ES) and Human Disturbance Intensity (HDI)

The LISA cluster map of bivariate local Moran describes the clustering patterns of ES and HDI (Figure 9). The clustering patterns of provisioning services and HDI are relatively scattered, and H_{(provisioning services)}H areas are widely distributed in the periphery of the urban center. The L_{(provisioning services)}H areas are most concentrated in the city center (Hong Kong, Macao, Shenzhen, and Guangzhou; Figure 9a). H_{(regulating and maintenance services)}H areas are mainly clustered in the coastal area, and in parts of Foshan and Huizhou. The L_{(regulating and maintenance services)}H areas are most concentrated in the center of the GBA (Figure 9b). For cultural services, the H_{(cultural services)}H areas are distributed across a large area and concentrated in the central cities of the GBA; L_{(cultural services)}H goes around the periphery of central cities (Figure 9c). There are similarities in the cluster patterns of LL and HL among the three ES groups. They are scattered across the edge of the GBA with different degrees of aggregation; the most obvious aggregation is for cultural services, followed by regulating and maintenance services, and provisioning services. The spatial patterns of insignificant areas are similar in all LISA cluster maps; they are mainly distributed in the periphery of the central cities due to relatively fewer human activities (e.g., population aggregation, land development, and infrastructure construction) in these areas in comparison to the metropolitan area of the GBA.

4. Discussion

4.1. Mechanism of Interactions between Human Disturbance and Ecosystem Services

A large number of job opportunities, platforms, and policies for supporting business startups and innovation, and complete infrastructural facilities for living and production attract people to the GBA. However, the GBA is facing the problem of increased differentiation of population quantity. Of the 1.5 million population increment in the Pearl River Delta region in 2018, Guangzhou and Shenzhen account for 0.40 million and 0.49 million, respectively. The populations of Foshan and Zhuhai are also growing notably, whereas the remaining five cities only saw population increments in the tens of thousands [30].

Most natural land has been converted to built-up land to meet the demand of urbanization and industrialization. For many years, coastal areas and offshore islands have been subjected to continuous reclamation projects, and terrestrial and marine resources have been utilized with high intensity [51]. Built-up land covers an area of 6769.17 km², which accounts for 12.18% of the total land area. It is mainly distributed over Guangzhou, Shenzhen, Hong Kong, and Macao.

The GBA has a complete transportation infrastructure and a dense road network in a radial distribution from the central cities to the surrounding areas. Foshan and Dongguan are making efforts to build a modern three-dimensional comprehensive transportation hub and accelerate the flow of people, logistics, capital, and information to the second-and third-tier cities; the first-tier cities are developing a hub-based economy in an effort to promote regional balance [52].

Shenzhen, Dongguan, Hong Kong, and Macao have highly developed economies and industries that require large energy resources to maintain daily operations. Compared with the light brightness value of other regions and highly populated bay areas, the degree of connectivity of the night-time lights of the GBA is significantly higher than that of the Tokyo Bay Area or the Yangtze River Delta urban agglomeration. The night-time light blocks are completely integrated, which indicates that the degree of urban integration in the GBA has reached a high level.

It can be seen from the aforementioned analysis that population and traffic pressure are the dominant factors of HDI in Guangzhou, Shenzhen, Dongguan, Foshan, Jiangmen, and Huizhou; there are greater material needs here and they are convenient for the delivery of provisioning services. The high- and medium-high-capacity areas of regulating and maintenance services are mainly distributed in the coastal areas of Guangzhou, Shenzhen, Zhuhai, Zhongshan, Zhaoqing, and Huizhou; these areas are mainly affected by the pressure from dense population and intense land-use. Cultural services are highly developed in Guangzhou, Shenzhen, Hong Kong, and Macao and the coastal zone of Zhuhai and Dongguan, where educational resources and recreational facilities are in great demand due to high-density population and highly developed economies. There is a complex feedback mechanism between human interference and ES. Human activities have positive effects on provisioning services and cultural services. Conversely, they can weaken the regulating and maintenance services, which will in turn impose restrictions on the delivery of provisioning services and cultural services. Areas with abundant ES are more conducive to carrying out various social and economic activities than areas with poor ecological functions.

4.2. Implications for Regional Planning

Implementing human activities without adversely affecting the sustainability of ecosystems, and minimizing the pressure on ES has been an essential concern for decision-makers worldwide. The current planning for developing functional zones attempts to alleviate the ecological environmental risks by directly delineating ecological protection areas in areas with a high concentration of ES, ecologically sensitive and fragile areas, and prohibiting or restricting economic and social activities. However, the study demonstrates that there is significant correlation between ES and HDI in the GBA, and that the existing planning overlooks the positive feedback between human activities and ES. Therefore, we should consider ES and human activities holistically and consider their interactions to minimize the negative impact on the natural environment and maximize ES value. Areas with high population density and high capacity for provisioning services and cultural services can be considered as ecological development zones for making full use of potential ecological values under the premise of priority protection. Areas with high capacity for regulating and maintenance services should be regarded as key ecological protection areas, and the strictest ecological protection measures should be adopted. In L_{(regulating and maintenance services)}H areas, the intensity of urban construction should be appropriately controlled, or the land-use types used should be adjusted to create ecological landscapes (such as artificial wetlands and urban green spaces) and alleviate the impact of human activities on the environment. In L_{(provisioning services, cultural services)}L areas, appropriate ecological restoration projects can be adopted due to lack of good natural conditions and abundant resources. Appropriate construction can be carried out to enhance their development vitality.

Urban planning and ecological landscape planning should be focused on H_{(regulating and maintenance services)}H areas. From the perspective of spatial correlation, most high-capacity areas of regulating and maintenance services are under high HDI. The natural habitats with important ecological value in coastal regions (i.e., mangrove wetlands, coral reef wetlands, and seagrass beds) are affected by intense human activities (i.e., aquaculture, fishing, foreign trade, shipping port transportation, and leisure tourism). For example, the Halophila beccarii Asch seagrass bed in Tangjiawan, Zhuhai, can lessen the impact of waves, consolidate the coastline, and protect seabed geology. At present, the seagrass beds are rapidly decreasing owing to high-intensity human activities [53]. National nature reserves such as Xiangtou Mountain in Huizhou and Daya Bay in the South China Sea have also become tourist attractions. A high HDI will cause the loss of regulating and maintenance services and lead to further weakening of other ES. Hence, it is necessary to coordinate the scope of human activities and minimize adverse effects on natural ecosystems. The intensity and types of human activities in areas with high regulating and maintenance services value should be strictly controlled to avoid degradation of ecological resilience.

Human disturbance not only affects ES in the area where disturbance occurs but it will also lead to changes of ecological value in surrounding areas through the exchange of energy and materials with neighboring regions via natural processes, such as hydrological and atmospheric movements [54,55]. The LISA cluster map demonstrates that H_{(provisioning services, cultural services)}H agglomerated in the center city of the GBA and large areas of L_{(provisioning services, cultural services)}H cluster in its the neighboring cities, which indicates that there is a negative spillover effect in the spatial relationship between ES and HDI [56]. The high-capacity region for provisioning services and cultural services can lose its functions due to surrounding intense human disturbance. Therefore, we should avoid environmentally destructive activities around key ecological protection areas. Coastal wetlands have high capacity for regulating and maintenance services, whereas coastal zones are also H_{(regulating and maintenance services)}H areas adjacent to the L_{(regulating and maintenance services)}H areas. Intense human activities inevitably affect ES directly or indirectly in ecologically valuable areas. Therefore, land development and economic activities should be carefully planned in its surrounding areas.

4.3. Advances and Limitations

This study fills the knowledge gap of the interaction between human disturbance and coastal ecosystems. Compared with the global-scale research in inland regions, our study is more specific and targeted. We fully considered the characteristics of coastal areas, such as coastal wetlands and navigable waterways when selecting indicators of human disturbance factors. ES provided by coastal terrestrial ecosystems, wetland ecosystems, and nearshore ecosystems were also comprehensively considered, and the kernel density of ES sources were used to visualize the spatial distribution of ES hotspots more objectively compared with the calculation of ecological service values. Geographically weighted Pearson correlation coefficient and the bivariate local Moran can spatially visualize and quantitatively reflect the correlation between ES and HDI in local space compared with the previous methods (like linear regression).

There are still some limitations. Many human activities affect ES; however, we evaluated the regional HDI through only four types of human pressure. More specific human activities (i.e., consumption habits and types of industry, tourism pressure, types of energy consumed, and other specific mechanisms) that exert a considerable influence on ecosystems should also be explored. Additionally, we used 3 × 3 km grids as spatial units to study the spatial relationship between ES and HDI; this may generalize or overlook some interactions that occur at smaller scales. Therefore, future research should use multiscale grid units or administrative units for evaluation in order to obtain more comprehensive and reliable results.

5. Conclusions

We quantitatively explored and visualized the spatial correlation between human disturbance intensity and ecosystem services in the Guangdong–Hong Kong–Macao Greater Bay Area. The major findings are summarized as follows.

The central and southern regions of the Guangdong–Hong Kong–Macao Greater Bay Area are under intense human disturbance. The kernel density of ecosystem service sources can spatially and objectively manifest the capacity degree of ecosystem services. Provisioning services mainly originate from the outer edge of the central cities. Cultural services are highly concentrated in the central cities of the study area. Regulating and maintenance services have higher density in the outermost peripheries of the Guangdong–Hong Kong–Macao Greater Bay Area. High-capacity regions of ecosystem services are mainly under high intensity of human disturbance. Regions with high provisioning service values are mainly under the pressures of population and traffic, whereas regulating and maintenance services and cultural services are mainly dominated by high-density populations; regulating and maintenance service values are also affected by land-use pressure.

The global correlation analysis shows that provisioning services are significantly negatively correlated with human disturbance intensity, whereas a positive correlation is strongest between cultural services and human disturbance intensity, followed by that between regulating and maintenance services and human disturbance intensity. However, the correlation between human disturbance intensity and ecosystem services is spatially nonstationary when examined locally. This spatial heterogeneity is influenced by specific human activities, natural conditions, and management mechanism in local area.

Human disturbance has negative spillover effects on ecosystem services, which requires avoiding destructive activities around key ecological protection areas. It is of great significance to identify the spatial correlation between ecosystem services and human disturbance intensity for maintaining ecological resilience, optimizing coastal ecological planning, and improving resources allocation to ensure more equitable access among areas.