Effects of Coastal Urbanization on Habitat Quality: A Case Study in Guangdong-Hong Kong-Macao Greater Bay Area

Abstract

Coastal areas are usually considered as pioneering areas for economic development and reform due to their unique geographical locations and ecological conditions. Correspondingly, rapid urbanization in coastal urban agglomerations has resulted in population concentration and land use/cover change (LUCC), leading to the decline of habitat quality and biodiversity. However, few studies have quantitatively explored the impacts of urban agglomeration expansion in coastal zones on habitat quality. Taking the Guangdong-Hong Kong-Macao-Great Bay Area (GBA) as a case study, we applied the Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model to assess the habitat quality based on land use data obtained from remote sensing images from 1980 to 2020 and developed a geographically weighted regression model to quantitatively analyze the response of habitat quality to urbanization. The results showed that (1) LUCC in the GBA was dramatic from 1980 to 2020, dominated by the shift from various land use types to construction land, which led to increased fragmentation; (2) the overall habitat quality showed a decreasing trend, characterized by low values in the central part and high values in the surrounding area; (3) population and construction land such urbanization elements had a more significantly negative effect on habitat quality changes, while the relationships among slope, road distance, and habitat quality changes were complex. Based on above analysis, this paper suggests that future land management in the GBA should develop in the direction of intensification, refinement, and regional integration.

1. Introduction

Coastal zones are interaction areas between terrestrial and marine ecosystems [1]. Their unique geographical locations and abundant natural resources make them important biological habitats, as well as key zones for urbanization and economic development [2]. With the continuous advancement of urbanization, land use changes caused by population concentration and industrial and infrastructure construction have led to habitat fragmentation and loss, thereby affecting biodiversity [3,4,5,6]. Furthermore, this trend has been getting worse over the past years [7]. Previous studies have suggested that coastal cities are expanding exponentially [8] and the loss of ecosystem services is particularly prominent in coastal zones [9]. Therefore, it is necessary to understand the state and variation process of habitat quality in coastal agglomerations and then to analyze the mechanism of habitat quality degradation due to urban expansion [10].

Land use and land cover change (LUCC) is the main threat to biodiversity loss, which is closely related to habitat quality decline [3,4,11]. Habitat quality is a core indicator of eco-environmental level [12]. Relevant studies are mainly based on LUCC, investigating the spatial and temporal evolution of habitat quality [13], taking scene simulations using urban expansion models (e.g., the Conversion of Land Use and its Effects at Small Region Extent model and the Cellular Automata model) [14,15], constructing regional ecological security patterns [16], and predicting the future direction of ecological patterns [17]. Researchers have found that the unpromising LUCC trend (i.e., the rapid increase of construction land and significant decrease of wetlands) has a close relationship with socio-economic development [1,18]. What deserves to be noticed is that this pattern is the result of rapid urban expansion [1]. In addition, urbanization also results in the decrease of landscape connectivity and the increase of fragmentation [19,20]. Studies on habitat landscape changes in cities are critical for biodiversity conservation [19,21]. Landscape analysis focuses on the type, shape, distribution, and arrangement of ecosystem components under different land cover conditions, which can reveal the urbanization effect on habitats inspatial dimension [14,18,22]. Different from other methods for monitoring variations in landscape patterns, i.e., land use models and remote sensing, the landscape index can provide precise information about landscape configuration and patterns [14,23]. It is a common method for quantitatively describing landscape fragmentation, diversity, and structure. However, while most researchers have explored the impacts of LUCC on habitat quality, study on spatial heterogeneity at a landscape level is scarce [22]; This issue needs to be addressed [22,24,25].

Habitat quality refers to the capacity of an ecosystem to provide the conditions for species to survive and multiply, reflecting the level of ecosystem services as an important aspect of biodiversity [26,27,28]. Habitat quality assessments have become an important component in assessments of ecosystem services [7]. They can help us understand the carrying capacity of regional ecosystems and the welfare gained by humans, making it possible to construct a sustainable management pattern [5,29]. Calderon and An analyzed the effects of pool, riffle, and run ecosystems on fish [30]; Cuffney et al. assessed the responses of benthic macroinvertebrates to urbanization and LUCC by field reconnaissance and according to various indicators [31]. Previous studies on biodiversity have mainly been conducted in the form of field surveys. Due to its high cost in terms of time and labor, this method is not suitable for large-scale studies [32]. In recent years, with progress in computer technology, GIS, and remote sensing technology, more attention is being paid to the dynamic monitoring of habitat quality on a large scale [12,25]. In this regard, there are two common methods: the first adopts an index system, e.g., the analytical hierarchy process, fuzzy comprehensive evaluation, and artificial neural networks [33,34,35]; the second includes ecosystem models, such as the Habitat Suitability Index model, Social Values for Ecosystem Services model, Integrated Valuation of Ecosystem Services and Trade-offs (InVEST) model [36,37,38], etc. Among them, the InVEST model generates habitat quality maps based on LUCC and the distribution of threat factors [39,40]. It is widely used to quantify habitat quality in different regions and on different scales due to its ease of use and flexible output [41,42].

Rapid coastal urbanization in recent years has profoundly affected land use patterns, resulting in the reduction of ecological land; this trend will continue [7]. Therefore, the influence of human activities, especially urbanization, on coastal ecosystems has become a hot spot in recent years. Numerous studies have explored the impact of coastal urbanization on regional ecosystems from the perspectives of ecological patterns, land use, and landscapes. Mayer-Pinto et al. investigated the structure and function of microbial communities in artificial and natural coastal environments [43]. They focused on how the differences between the two ecosystems were reflected in biodiversity changes. Dupras et al. studied the variation of ecosystem services on the Mediterranean coast based on LUCC [44]. Aguilera et al. investigated the effects of urbanization on coastal ecosystem spatial connectivity from a landscape pattern perspective [45]. More and more studies have shown that urbanization has led to coastal ecosystem degradation. Therefore, land use optimization seems to be necessary in urban construction. Common methods with which to analyze the influence of urbanization on habitat quality include correlation analysis [46], regression model [45], and index appraisal [1]. These approaches are appropriate for solving questions about the impact of a single factor (e.g., land use conversion) on habitat quality but are not suitable for analyses of the extent and direction of natural and social factors in the process of urbanization. Some studies have noted that habitat loss caused by LUCC is the result of interactions between natural ecosystems and humans. Such influences of LUCC are usually spatially autocorrelated because of their similar geographical environments [47]. The geographically weighted regression (GWR) model introduces the location attributes into regression analysis; it is a robust method for solving problems about spatial heterogeneity [48,49]. Currently, the GWR model is being used in various spatial correlation studies. Su et al. analyzed the spatial changes between agriculture landscape and urbanization; those authors noted that the GWR model has better explanatory accuracy and yields richer information [50]. Shearmur et al. also used it to explore the spatial and temporal differences behind the reasons for employment growth in Canada, showing that the model has strong explanatory power for socioeconomic problems [51].

The Guangdong-Hong Kong-Macao Greater Bay Area (GBA), which borders on the South China Sea, is one of the most economically open and active regions in China. Over the past four decades, it has experienced rapid urbanization and dramatic LUCC, which has placed increasing pressure on regional ecological areas and biodiversity [52,53,54]. To this end, taking GBA as a case study, the aim of this study is to fill in the gap in habitat quality assessments during the period of rapid coastal urbanization. The research objectives of this study are: (1) to analyze the rate and direction of land use changes in GBA; (2) to investigate the landscape pattern features and variation tendency; (3) to explore the spatial-temporal evolution characteristics of habitat quality from 1980 to 2020; and (4) to quantitatively reveal the impact of coastal urbanization on habitat quality based on the dimensions of nature and society.

2. Materials and Methods

2.1. Study Area

The GBA (21°17′ N–23°55′ N, 111°59′ E–115°25′ E) is a cluster of nine cities in Guangdong Province, namely, Guangzhou, Jiangmen, Zhaoqing, Zhuhai, Shenzhen, Huizhou, Dongguan, Foshan, and Huizhou, as well as the two special administrative regions of Hong Kong and Macao. It is located in southern China, at the intersection of the Pearl River and the South China Sea, with a long coastline and vast maritime area (Figure 1). The GBA covers approximately 56,000 km². The east, north, and west sides are surrounded by hills, with plains in the center and southeast. GBA is the fourth largest bay area in the world after the New York Bay Area and the San Francisco Bay Area in USA and the Tokyo Bay Area in Japan. As such, it is an important carrier for the country to construct a world-class city cluster and participate in global enterprise. In 2020, the total population of the GBA exceeded 86 million, the regional gross domestic product (GDP) reached 166.88 billion dollars, and the urbanization rate reached over 80% (Guangdong Provincial Bureau Statistics and Survey office of the National Bureau of Statistics in Guangdong, 2021). The “Outline Development Plan for the GBA”, released in 2019, pointed out the need to develop an international first-class bay area for living, working, and travelling, which underlined the urgency of the eco-environmental protection.

2.2. Data Sources and Preparation

The datasets in this study mainly included the following three aspects:

• Land use data: Land use data of study area for 1980, 1990, 2000, 2010, and 2020 were obtained from the Land use and Land Cover of China (CNLUCC) data provided by the Resource Environment Science and Data Center of the Chinese Academy of Sciences (https://www.resdc.cn/, accessed on 1 July 2022), with a spatial resolution of 30 m. CNLUCC data based on Landsat TM series remote sensing images were processed using the supervised classification method, with an accuracy of more than 90% [55]. According to the present research needs and the Chinese land resource classification system, we reclassified the original land cover types into six categories: cultivated land, forestland, grassland, water bodies, construction land, and unused land. The spatial distributions of land use types in the GBA from 1980 to 2020 are displayed in Figure 2.
• Socioeconomic data: Raster datasets of population density in 1990, 2000, and 2010 were obtained from the same source as above, with a spatial resolution of 1 km (Figure 3a–c). For the missing data, we calculated the weights of each district and county, referred to statistical yearbooks, and completed them according to adjacent years. The vector data of administrative boundaries and roads came from the National Geomatics Center of China (https://nfgis.nsdi.gov.cn/, accessed on 1 July 2022). Considering that the changes of roads were not obvious during a certain period, we only chose Road I and II in two periods (Figure 3d), corresponding to the periods of 1980–2000 and 2000–2020, respectively.
• DEM data: SRTM elevation data, with a spatial resolution of 30 m, were obtained from the official website (https://earthexplorer.usgs.gov/, accessed on 1 July 2022). The slope data of study area were generated from DEM, as shown in Figure 3e.

2.3. Methods

2.3.1. Analysis of Land Use Change

Land use transfer matrices can quantitatively describe the land use type and status in a certain period. We used the overlay analysis tool in ArcGIS to calculate the land use transfer matrix for different periods and analyzed the intensity of land use changes based on the land use dynamic degree.

A single land use dynamic degree represents an increase or decrease of a specific land use type area during the study period [56]. The calculation formula is as follows:

$$\mathrm{K}=\frac{{\mathrm{LU}}_{\mathrm{b}}-{\mathrm{LU}}_{\mathrm{a}}}{{\mathrm{LU}}_{\mathrm{a}}}$$

(1)

where ${\mathrm{LU}}_{\mathrm{a}}$ and ${\mathrm{LU}}_{\mathrm{b}}$ denote an area of the same land use type at the beginning and end of the study period, respectively.

Comprehensive land use dynamic degree (LC) quantitatively describes the land use change rate in a particular period, reflecting comprehensive changes in land use area [32,57]. The formula is as follows:

$$\mathrm{LC}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}∆{\mathrm{LU}}_{\mathrm{ij}}}{2{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{LU}}_{\mathrm{i}}}\times \frac{1}{\mathrm{T}}\times 100\%$$

(2)

where ${\mathrm{LU}}_{\mathrm{i}}$ represents the land use area at the initial time and $∆{\mathrm{LU}}_{\mathrm{ij}}$ represents the land use area converted from type $\mathrm{i}$ to type $\mathrm{j}$ within period $\mathrm{T}$.

2.3.2. Landscape Index

It is essential to consider whether the indices are independent when selecting them. In order to describe the fragmentation, aggregation, dominance, and diversity of landscapes, we selected indices at the class and landscape levels, as shown in

Table 1. Landscape indices used in this study.

	Landscape Metrics	Formula	Description
Class level	Number of patches (NP)	$\mathrm{NP}=n_i$, where ${\mathrm{n}}_{\mathrm{i}}$ is the number of patches of type i	$\mathrm{NP}\in [1,+\infty )$ reflects the spatial pattern of a landscape and has a positive correlation with fragmentation.
	Mean patch size (MPS)	$\mathrm{MPS}=\frac{{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{ij}}}{{\mathrm{n}}_{\mathrm{i}}}$, where ${\mathrm{a}}_{\mathrm{ij}}$ is the area of patches of type i	$\mathrm{MPS}\in (0,+\infty )$; the lower the value, the greater fragmentation of patch class.
	Largest path index (LPI)	$\mathrm{LPI}=\frac{\max \left({\mathrm{a}}_{\mathrm{ij}}\right)}{\mathrm{A}}\times 100\%$, where $\mathrm{A}$ is total area	$\mathrm{LPI}\in (0,100]$ refers to the ratio of the largest patch area to the total landscape area. Its value can characterize the dominance of a given landscape type.
Landscape level	Patch density (PD)	$\mathrm{PD}=\frac{\mathrm{N}}{\mathrm{A}}$, where $N$ is the is the total number of patches in the landscape	$\mathrm{PD}\in (0,+\infty ]$; a higher value means more patches per unit area and greater fragmentation.
	Perimeter-area fractal dimension (PAFRAC)	$\mathrm{PAFRAC}=\frac{2}{\frac{\mathrm{N}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}\left(\ln {\mathrm{p}}_{\mathrm{ij}}\times \ln {\mathrm{a}}_{\mathrm{ij}}\right)-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{ij}}\times {{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{a}}_{\mathrm{ij}}}{\mathrm{N}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{ij}}{}^2-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{n}}\ln {\mathrm{p}}_{\mathrm{ij}}}}$	$\mathrm{PAFRAC}\in [1,2]$ describes the complexity of the landscape. When the value tends to 1, it means that the patch shape is simple.
	Aggregation index (AI)	$\mathrm{AI}=\frac{{\mathrm{g}}_{\mathrm{ii}}}{\max \to {\mathrm{g}}_{\mathrm{ii}}}$, where ${\mathrm{g}}_{\mathrm{ii}}$ is the number of similar neighboring patches of a given landscape type	$\mathrm{AI}\in (0,100]$ examines the degree of clustering of classes within a landscape; the smaller the value, the higher the dispersion degree.
	Shannon’s diversity index (SHDI)	$\mathrm{SHDI}=-{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}{\mathrm{P}}_{\mathrm{i}}\ln {\mathrm{P}}_{\mathrm{i}}$, where ${\mathrm{P}}_{\mathrm{i}}$ is the probability of the occurrence of landscape patch type i in the landscape	$\mathrm{SHDI}\in [0,+\infty )$ reflects the richness of each patch type within the landscape.

[58], which were calculated in the Fragstats 4.2 software.

2.3.3. Habitat Quality Assessment Model

The InVEST habitat quality model assesses habitat quality and the threat degree to biodiversity based on LUCC, which is reflected in the habitat quality index, ranging from 0 to 1; the larger the value, the better the habitat quality. The formula is as follows:

$${\mathrm{Q}}_{\mathrm{xj}}={\mathrm{H}}_{\mathrm{j}}\left[1-\left(\frac{{\mathrm{D}}_{\mathrm{xj}}^{\mathrm{z}}}{{\mathrm{D}}_{\mathrm{xj}}^{\mathrm{z}}+{\mathrm{k}}^{\mathrm{z}}}\right)\right]$$

(3)

where ${\mathrm{Q}}_{\mathrm{xj}}$ is the habitat quality of grid $\mathrm{x}$ in land use type $\mathrm{j}$; ${\mathrm{H}}_{\mathrm{j}}$ is the habitat adaptability of land use type $\mathrm{j}$; $\mathrm{z}$ and $\mathrm{k}$ are the normalization constant and half-saturation constant, respectively; and ${\mathrm{D}}_{\mathrm{xj}}$ is the habitat degradation degree, defined as follows:

$${\mathrm{D}}_{\mathrm{xj}}={\displaystyle \sum }_{\mathrm{r}=1}^{\mathrm{R}}{\displaystyle \sum }_{\mathrm{y}=1}^{{\mathrm{Y}}_{\mathrm{r}}}\left(\frac{{\mathsf{\omega }}_{\mathrm{r}}}{{{\displaystyle \sum }}_{\mathrm{r}=1}^{\mathrm{R}}{\mathsf{\omega }}_{\mathrm{r}}}\right){\mathrm{r}}_{\mathrm{y}}{\mathrm{i}}_{\mathrm{rxy}}{\mathsf{\beta }}_{\mathrm{x}}{\mathrm{S}}_{\mathrm{jr}}$$

(4)

where $\mathrm{R}$ and ${\mathrm{Y}}_{\mathrm{r}}$ indicate the number of threat factors and grids, respectively; ${\mathsf{\omega }}_{\mathrm{r}}$ is the weight of threat factor $\mathrm{r}$; ${\mathrm{r}}_{\mathrm{y}}$ is the threat value of grid $\mathrm{y}$; ${\mathrm{i}}_{\mathrm{rxy}}$ indicates the impact of threat value ${\mathrm{r}}_{\mathrm{y}}$ in grid $\mathrm{y}$ on habitat grid $\mathrm{x}$; ${\mathsf{\beta }}_{\mathrm{x}}$ is the anti-interference level of the habitat; and ${\mathrm{S}}_{\mathrm{jr}}$ is the relative sensitivity of habitat type $\mathrm{j}$ to threat factor $\mathrm{r}$.

There are two attenuation methods for ${\mathrm{i}}_{\mathrm{rxy}}$, represented by the following equations:

$$\mathrm{Linear} \mathrm{attenuation}:{ \mathrm{i}}_{\mathrm{rxy}}=1-\left(\frac{{\mathrm{d}}_{\mathrm{xy}}}{{\mathrm{d}}_{\mathrm{r} \max }}\right)$$

$$\mathrm{Exponential} \mathrm{decay}:{ \mathrm{i}}_{\mathrm{rxy}}=\exp \left[-\left(\frac{2.99}{{\mathrm{d}}_{\mathrm{r} \max }}\right){\mathrm{d}}_{\mathrm{xy}}\right]$$

where ${\mathrm{d}}_{\mathrm{xy}}$ is the distance between habitat grid $\mathrm{x}$ and threat grid $\mathrm{y}$ and ${\mathrm{d}}_{\mathrm{r} \max }$ is the maximum impact distance of threat factor $\mathrm{r}$.

The input parameters of the InVEST model include land use maps for each period, threat factor data (including maximum impact distance, weight, and degradation type), threat factor maps, sensitivity of land use types to each threat, and the level of legal protection from disturbance in each grid. In this study, we considered all grids in the study area to be equally protected and mainly considered the effects of other parameters.

With this model, it is first necessary to distinguish the habitats. As reflections of human activity, construction land and cultivated land can both lead to habitat fragmentation, threatening the natural ecosystem to a certain extent. Therefore, we defined cultivated land, construction land, and unused land with low vegetation coverage such as sand, bare ground, and bare rocky land as threat sources, and forestland, grassland, and water bodies as habitats.

The parameters to be adjusted according to the specific conditions in the study area include the maximum impact distance (i.e., the weights) of threat factors, the suitability scores of different land use types, and the sensitivity to threat sources, which were determined based on the model guidebook and research results from other scholars [16,39,59,60]. The specific parameters are shown in the

Table 2. Threats and their maximum influence distances and weights.

Threat Factors	${\mathit{d}}_{\mathit{r} \mathit{m}\mathit{a}\mathit{x}}\mathbf{/}\mathbf{km}$	Weight ωr	Distance-Decay Function
CT	8	0.5	linear
UL	10	1	exponential
RS	8	0.8	exponential
OCL	9	0.9	linear
UL	5	0.3	exponential

and

Table 3. Habitat suitability and sensitivity of each land use type.

Land Use Types	Suitability	CT	UL	RS	OCL	UL
CT	0.5	0	0.7	0.5	0.6	0.3
WL	1	0.7	0.9	0.7	0.8	0.3
SH	1	0.6	0.8	0.6	0.7	0.2
SP	1	0.5	0.7	0.5	0.6	0.2
OWL	0.8	0.4	0.7	0.5	0.6	0.2
GL	0.9	0.5	0.8	0.7	0.8	0.4
RV	1	0.7	0.9	0.7	0.8	0.2
LK	1	0.7	0.9	07	0.8	0.2
RE	1	0.3	0.5	0.3	0.4	0.2
TD	1	0.8	1	0.8	0.7	0.2
FA	0.8	0.8	1	0.8	0.7	0.3
UL	0	0	0	0	0	0
RS	0	0	0	0	0	0
OCL	0	0	0	0	0	0
UL	0.3	0.3	0.5	0.4	0.5	0
SA	1	0.2	0.6	0.4	0.5	0.1

CT: cultivated land; WL: woodland; SH: shrubbery; SP: sparse woodland; OWL: other woodland; GL: grassland; RV: rivers; LK: lakes; RE: reservoir; TD: tidal; FA: flat; UL: Urban land; RS: rural settlement; OCL: other construction land; UL: unused land; SA: sea.

.

2.3.4. Quantitative Analysis of the Impact of Urbanization

1.

Fundamental principle

Regression analysis is a common method to explore relationships among variables. As a global regression, ordinary least squares (OLS) is based on the minimum estimated residual sum of squares to determine the regression model. It can be stated as:

$$\mathrm{y}={\mathsf{\beta }}_0+{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{k}}{\mathsf{\beta }}_{\mathrm{i}}{\mathrm{x}}_{\mathrm{i}}+\mathsf{\epsilon }$$

where $\mathrm{y}$ is the explanatory variable; ${\mathsf{\beta }}_0$ is the intercept; ${\mathsf{\beta }}_{\mathrm{i}}$ is the regression parameter coefficient of independent variable ${\mathrm{x}}_{\mathrm{i}}$; $\mathsf{\epsilon }$ is the random error term fitted to the normal distribution; and $\mathrm{k}$ is the number of independent variables.

Due to the different natural geographic conditions, regulatory policies, and effects across the study area, the impacts of urbanization on habitat quality are complex and changeable. In other words, the drivers of urbanization are spatially heterogeneous. The GWR model introduces the location attributes as distances based on weights into the regression equation, which can describe the spatial non-stationary nature of urbanization factors more precisely and make the results more reliable. Therefore, with each regression point as a center, this study established local regression equations for GWR analysis within a certain radius. It can be expressed as:

$${\mathrm{y}}_{\mathrm{i}}={\mathsf{\beta }}_0\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)+{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{k}}{\mathsf{\beta }}_{\mathrm{j}}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right){\mathrm{x}}_{\mathrm{ij}}+{\mathsf{\epsilon }}_{\mathrm{i}}$$

where $\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ is the geographic coordination for location $\mathrm{i}$; ${\mathsf{\beta }}_{\mathrm{j}}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ represents the regression coefficient for independent variable ${\mathrm{x}}_{\mathrm{j}}$ at location $\mathrm{i}$; ${\mathsf{\beta }}_0\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ is the intercept; and ${\mathsf{\epsilon }}_{\mathrm{i}}$ is the random error term where ${\mathsf{\epsilon }}_{\mathrm{i}}\sim \mathrm{N}\left(0,{\mathsf{\sigma }}^2\right)$.

The regression coefficient at position $\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ is given by $\hat{\mathsf{\beta }}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$:

$$\hat{\mathsf{\beta }}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)={\left({\mathrm{X}}^{\mathrm{T}}\mathrm{W}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)\mathrm{X}\right)}^{-1}{\mathrm{X}}^{\mathrm{T}}\mathrm{W}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)\mathrm{Y}$$

where $\mathrm{W}\left({\mathrm{u}}_{\mathrm{i}},{\mathrm{v}}_{\mathrm{i}}\right)$ is the spatial weight matrix.

The calculation of weights depends on the distance decay function. Gaussian function is one of the commonly used methods, which can be stated as:

$${\mathsf{\omega }}_{\mathrm{ij}}=\exp {\left(-\frac{{\mathrm{d}}_{\mathrm{ij}}}{\mathrm{b}}\right)}^2$$

where ${\mathrm{d}}_{\mathrm{ij}}$ is the distance between sample points $\mathrm{i}$ and $\mathrm{j}$; and $\mathrm{b}$ is the bandwidth.

Bandwidth is the distance band for each local regression equation. It directly affects the fitting accuracy of the model and is perhaps the most important parameter to consider for GWR [61]. This study adopted the AICc method to determine the optimal bandwidth with the minimum AIC value.

2.

Model construction

Natural elements and human activities are key factors contributing to habitat quality differentiation [53]. Firstly, terrain is generally considered to affect the density, scale, and spatial distribution of urban land [62]. Since the terrain of the study area showed obvious ring-shaped features, we considered slope as a geographical determinant of urban expansion when analyzing the direct influence of natural factors on habitat quality. Secondly, socioeconomic indicators reflect the intensity of human activities, which indirectly affect habitat quality. Among them, land urbanization, population urbanization, and transportation are the driving force, core, and carrier of urbanization, respectively [63]. Therefore, we selected the construction land area, population, and road distances as socioeconomic factors. In other words, taking the changing values of habitat quality as dependent variables, the slope, changes in construction land area, population, and road distances in different periods were chosen as explanatory variables to establish the regression model.

In order to unify the resolution, we used ArcGIS to create fishing grids with a spatial resolution of 1 km as basic units for the statistics of urban drivers and habitat quality changes in the GBA. “Near Analysis” and “Zonal Statistics” were used to calculate the average slope in each grid, changes of construction land area and population in different periods, the nearest distance from the center of a grid to an arterial road, as well as the changes of average habitat quality value in each grid. Finally, this raster layer, with a resolution of 1 km, was used as the input of the regression model. Due to the lack of data in Hong Kong and Macao, this part only considered the nine cities in the Pearl River Delta. Meanwhile, in order to reduce the influence of errors in the supplementary data, we compared the OLS and GWR results for the 1980–2000 and 2000–2020 periods to choose the optimal model.

We checked whether there were any residual explanatory variables according to the OLS output report tables generated by ArcGIS. The results showed that the VIF value of each variable was less than 10 and passed the significance test; therefore, all parameters could be used in the GWR model. In this case, we used the geographically weighted regression tools in ArcGIS and chose the AICc bandwidth method. The fitting results are shown in

Table 4. Comparison of two regression models.

Measuring Metrics	1980–2000		2000–2020
	OLS	GWR	OLS	GWR
Sigma	0.057	0.050	0.072	0.056
AICc	−152,720.269	−167,632.616	−128,410.083	−154,192.882
R2	0.683	0.765	0.778	0.867
Radj2	0.683	0.762	0.778	0.865

and Figure 4.

It can be considered that the GWR model is a significant improvement on the OLS model when the difference of AICc values is greater than 3. The above table shows that GWR model has lower Sigma and AIC values and better goodness-of-fit, so the regression coefficient results of the GWR model were used for the impact analysis. The local R² results for this model are shown in Figure 5. It can be seen that compared with the period of 1980–2000, the range of better fitting results from 2000 to 2020 is wider, especially in the northwest and southwest margins of the GBA. This is mainly because the changes in habitat quality were not significant in these regions from 1980 to 2000 (as discussed in detail in Section 3.3 and Section 3.4). However, in general, the selected variables are appropriate for constructing the GWR model.

3. Results

3.1. Land Use Change Characteristics

We calculated the land use dynamics and the transfer matrix in the four periods of 1980–1900 (I), 1990–2000 (II), 2000–2010 (III), and 2010–2020 (IV). The results are shown in Figure 4.

It can be seen that cultivated land, forestland, and construction land have always been the main land types in the GBA, accounting for about 90% of the total area (Figure 6a). As shown, the area proportion of each land type has undergone great changes over the years. Notably, the area of construction land has increased by 5531.29 km², with the ratio of 4.73% in 1980 increasing to 14.79% in 2020; this land types is mainly distributed in the central and southeastern area. Its growth accelerated annually, reached the peak of 6.01% from 2000 to 2010, and then slowed down. Correspondingly, the cultivated land area showed on a general downward trend, but the rate was controlled in the last period. Although the area of unused land is decreasing at a relatively high speed year by year, the change is not significant in the whole study area due to the low area ratio. The water area was reduced after 2000 but increased slightly over the whole period. The areas of forestland and grassland were relatively stable. LC was on the rise from 1980 to 2010, increasing by 1.65 times (Figure 6b). This indicates that the rate of land use in the GBA increased significantly in this period and then slowed down. In addition, after 1990, the LC was greater than the overall dynamic degree during the study period (0.21%).

The area of cultivated land declined by 1186.74 km² from 1980 to 1990; it was mainly converted into construction land, water bodies, and forestland (

Table 5. Land use transfer matrix by period.

Period	Land Use Type	CT	FL	GL	WB	CL	UL
I (/km2)	CT	15,424.87	292.76	16.25	401.83	473.60	2.29
	FL	235.33	30,333.68	38.57	42.84	81.50	0.24
	GL	15.30	171.13	1211.75	4.42	6.54	0.04
	WB	120.37	41.31	4.51	3164.17	26.53	0.10
	CL	76.16	30.94	2.37	18.68	2466.51	0.02
	UL	45.98	6.61	0.21	77.86	11.27	12.14
II (/km2)	CT	13,808.67	288.80	16.40	796.60	1015.63	0.41
	FL	285.58	30,156.38	49.06	53.60	349.69	0.31
	GL	17.75	61.31	1151.30	4.76	38.82	0.11
	WB	153.71	51.91	3.90	3387.48	151.77	0.21
	CL	76.28	36.84	2.77	16.01	2934.92	0.01
	UL	0.25	0.62	0.09	0.12	0.01	14.54
III (/km2)	CT	11,516.64	354.63	19.51	684.89	1766.25	0.27
	FL	307.93	29,378.12	50.87	95.44	762.81	0.60
	GL	21.32	102.73	1019.86	14.06	65.53	0.05
	WB	576.60	68.89	7.34	3090.56	516.04	0.12
	CL	183.08	135.39	6.34	79.46	4106.61	0.07
	UL	2.33	1.02	0.10	0.42	4.44	7.22
IV (/km2)	CT	11,404.05	203.63	21.65	186.64	790.19	0.14
	FL	226.06	29,156.00	86.05	104.85	456.87	0.29
	GL	13.74	36.46	1011.18	7.39	33.87	0.03
	WB	102.37	67.58	12.33	3509.94	270.58	0.14
	CL	337.45	180.03	47.48	82.57	6572.97	0.02
	UL	0.42	0.30	0.05	0.18	1.49	6.17

). From 1990 to 2000, the main trends were the transformation of cultivated land, forestland, and water bodies into construction land; this mainly occurred on the banks of the Pearl River (e.g., Foshan, Shenzhen, and Dongguan) (Figure 7). The scale of LUCC enlarged in this period. From 2000 to 2010, as the transformation of other land use types, especially cultivated land to construction land, continued and the level became more intense, the expansion rate of construction land reached a peak, displaying a tendency of expansion from the center to the periphery. Meanwhile, 576.60 km² of water bodies were converted into cultivated land. This mainly occurred in Foshan, accounting for 73.68% of the total area transferred in this way.

The Sankey diagram below clearly reflects the main direction and scale of land use conversion over the years (Figure 8). Notably, over the past 40 years, the increase in construction land mainly came from the occupation of cultivated land, forestland, and grassland. The expansion mainly occurred in coastal cities and their main districts such as Dongguan, Shenzhen, the Yuexiu District of Guangzhou, the Chancheng District of Foshan, etc. At the same time, there were also conversions among cultivated land, water bodies, and forestland, due to adjustments of the agricultural structure and the implementation of a farming compensation policy.

3.2. Land Use Change from the Perspective of Landscape Indices

3.2.1. Class Level Analysis

Figure 9 shows the results of the landscape pattern indices at the class level. NP and MPS indicate the fragmentation of landscape classes. NP of construction land showed a downward trend, even if the values were much higher than in other classes. Additionally, MPS was relatively small but obviously increased from 2010 to 2020. This trend indicates the continuous expansion and integration of construction land, developing in spatial congregation. From 1980 to 2010, cultivated land NP increased by 41.81%, and its MPS had the greatest change, with a decrease of 46.48%. NP and MPS of water, grassland, and unused land remained steady throughout the study period, although there were some fluctuations. LPI can indirectly reflect the level of human disturbance. Although forestland LPI fluctuated very little, the changes were notably more substantial than in other classes. This indicated that forestland was the dominant landscape at the class level, and that it was not significantly influenced by human activities.

3.2.2. Landscape Level Analysis

PD reflects landscape heterogeneity in a unit area. From 1980 to 2010, PD consistently increased from 0.459 to 0.503 (

Table 6. Landscape indices at the landscape level.

Indices	1980	1990	2000	2010	2020
PD	0.459	0.464	0.473	0.503	0.444
PAFRAC	1.447	1.449	1.450	1.407	1.408
SHDI	1.112	1.117	1.167	1.204	1.220
AI	96.684	96.626	96.508	96.576	96.731

). It was worth mentioning that PD decreased to 0.444 in 2020. The possible reason for this was that construction land gradually integrated and congregated from the scattered patches into large patches in that period. PAFRAC showed a peak of 1.450 in 2000, indicating the most complex patch shape at this time, before decreasing to 1.407 in 2010 and varying slightly in the other periods. SHDI is sensitive to the unbalanced distribution of landscape classes. From 1980 to 2020, SHDI increased steadily from 1.112 to 1.220, reflecting that the landscape class showed balanced development. AI decreased first then increased during the whole period. At the beginning of economic development, large-scale land resources were divided, leading to a scattered distribution. The combination and annexation of construction land patches was the main reason for the slight increase of AI after 2010.

3.3. Characteristics of Habitat Quality Changes in the Context of Urbanization

In reference to the results from and reclassification methods used in previous studies, we classified habitat quality values as poor (I: 0–0.3), medium (II: 0.3–0.6), good (III: 0.6–0.8), and excellent (IV: 0.8–1) and displayed them hierarchically (Figure 10).

Spatially, it showed that the habitat quality was lower in the middle of GBA and higher in its surroundings. The level of habitat quality differed significantly within the different cities. As the leading cities, the average values of the habitat quality index in Guangdong and Shenzhen were only 0.68 and 0.60, respectively. The habitat quality of cities along the Pearl River, including Dongguan, Foshan, Zhongshan, and Zhuhai, were in the medium level (Figure 11a). Peripheral cities in GBA, including Zhaoqing, Jiangmen, and Huizhou, had better habitat quality. It is worth noting that Hong Kong still had a high level, at 0.75 or even higher, because the city has been at complete urbanization phase.

Figure 9b reveals that the average habitat quality presented a downward trend from 1980 to 2020 in the GBA. The decrease speed accelerated over the years in the first three periods, and the area proportion with poor habitat quality increased from 5.12% to 14.96% (Figure 11c). In the areas with decreased habitat quality, 5.43% of them dropped by more than 0.5, while 94.57% dropped by 0–0.5. Four cities, i.e., Dongguan, Shenzhen, Foshan and Zhongshan, experienced the most significant decline, with decreases of 40.81%, 37.21%, 26.80%, and 19.93%, respectively.

3.4. Habitat Quality Response to Urbanization Factors

3.4.1. Impact of Natural Factors

The regression coefficients between slope and habitat quality differed significantly within the two phases (Figure 12). Generally, the extreme values of coefficients were greater than the marginal zones in the central and southeastern areas with less slope than the marginal zones. Initially, the significant negative correlation was mainly distributed in the central and southern cities, including Foshan, Zhongshan, and Dongguan. After 2000, the dependence of the slope coefficient on habitat quality changed to a positive correlation, and the range of significant negative impact expanded to the surrounding areas.

3.4.2. Impact of Socioeconomic Factors

The expansion of construction land mainly came at the expense of cultivated land, which was the chief reason for habitat quality reduction. The variation of habitat quality was significantly negatively correlated with construction land, which was distributed widely (Figure 13a). Although the previous regression result exhibited greater extreme coefficient values, the negative area was extended to the entire GBA from 2000 to 2020. Notably, most areas, and in particular, Shenzhen, Guangzhou, and Zhongshan, experienced significant habitat quality degradation due to the increasing impacts of construction. Population also had a negative effect on habitat quality. (Figure 13b) The negative correlation was most prominent in the cities with relatively low population densities and extensive vegetation coverage, such as Zhaoqing, Jiangmen, and Huizhou, because the habitat quality there was more sensitive to population growth. Transportation is a critical foundation for urbanization. Road construction usually results in the degradation and fragmentation of the ecological environment. The relationships between road distance and habitat quality were varied in different places (Figure 13c). From 1980 to 2000, the negative impact was concentrated in the cities with denser road webs such as Foshan, Zhongshan, and Shenzhen. Along with the improvement of transportation infrastructure, the road network became denser and denser, and the area with negative impact expanded further. This influence was more notable in Guangzhou, which has become an important transportation city.

4. Discussion

4.1. Spatial Pattern of Habitat Quality

The terrain conditions of the GBA shape the distribution pattern of habitat quality as “higher in the middle and lower in the surroundings”. In the vast central plains, frequent human activities and large-scale land reclamations have accelerated the disturbance of habitats, resulting in a lower level of habitat quality. Relatively complete forest and grassland resources are abundant in the western and northeastern areas with more mountains and hills, which provide a good basis for biodiversity. Moreover, the scale of land use change was small within the whole city; thus, the level of habitat quality was higher in the marginal cities than that in the central cities. Additionally, over a long period of time, cities of the east coast of the Pearl River, represented by Shenzhen, Dongguan, and Huizhou, developed faster than those on the west coast [64]. More fragmented patches with low habitat quality were concentrated in the southeastern hinterland. In a word, in areas with higher development levels and economic vitality, the negative impacts of human activities on the natural environment are greater, and thus, the decline in habitat quality is more significant. The same applies within cities. The most intensive socioeconomic activities are concentrated in central urban areas, and the patches with low habitat quality tend to spread around these areas.

4.2. Drivers of Habitat Quality Change

Our results are consistent with previous research, suggesting that during rapid urbanization, excessive human activities and the increase of demand for construction land to support development accelerate the evolution of land use patterns, directly resulting in significant habitat quality degradation [40,65,66]. In terms of indirect factors, based on previous studies, we consider that geographical conditions form the foundation of urbanization, as they determine the orientation of urban development to a certain extent [53,57]. Policy guides the development of urban agglomeration [67]. Population growth and road construction are the driving forces for urban expansion. In 1980, the pattern of opening up this coastal area, marked by the foundation of the Special Economic Zone, was initiated, and the road network in China began to expand from that time. Economic development and transportation construction mutually promoted each other, providing the conditions for population flow [68]. The massive influx of immigrants and the growth of the non-agricultural population were often associated with the rapid expansion of residential and commercial land. At this time, the region was still in the initial stage of economic construction. The rate of land use change grew with the increase in land demand. With the growth of industry and commerce and the stimulating effect of the real estate craze, various activities divided and occupied large swathes of cultivated land and forestland. The expansion of construction land came at the cost of land use types with higher habitat suitability. This process has destroyed the relative integrity and connectivity of land resources and led to habitat degradation. The period of 2000–2010 saw the most rapid urban development and economic growth. Excluding Hong Kong and Macao, the gross regional product of the nine cities in the Pearl River Delta increased from 847.129 billion yuan to 380.2846 billion yuan, and the population increased from 42.8978 million to 56.2295 million. The extent and speed of land use change reached a peak during this period. On the one hand, the agglomeration effect of urban development caused construction land to expand from the main urban area to the surrounding area, showing the trend of spatial concentration. On the other hand, non-agricultural construction and agricultural structure adjustments accelerated the fragmentation of cultivated land and forest patches. Forestland reclamation and the replacement of dike ponds with highly economically beneficial construction land were common. At the same time, in order to compensate for the cultivated land that was lost during the early industrial development, large-scale land reclamation from beaches or the sea led to a significant reduction of water bodies. This was considered one of the main causes of habitat quality decline in this period. Additionally, road construction has contributed to the population concentration, leading to landscape fragmentation and the reduction of biological habitats, which have further deepened the negative impact on habitat quality [14]. After 2010, economic growth tended to be stable and urban development reached a level of maturity. With an increase of governmental supervision and ecological regulations, urban planning has become more rational. Additionally, regional ecological degradation and habitat quality decline have been effectively limited. Especially since the implementation of the Pearl River Delta Reform and Development Plan, cities have intensified their land management and construction activities. The disordered sprawl of construction land has also been controlled. For Hong Kong and Macao, conservation areas, such as urban catchments and country parks, were delimited at the beginning of planning. Thus those cities still show a relatively stable level of habitat quality for high-density cities.

Our study provides a new perspective for analyzing the relationship between the urbanization process and changes in habitat quality caused by LUCC. Compared with previous studies, this paper takes spatial heterogeneity into consideration. Combined with a dynamic analysis of regional landscapes, it reveals the spatial characteristics and underlying process of urban expansion and clearly reflects the spatial variability of such impacts [22,47]. However, there are some limitations in this study. For example, to ensure the consistency of the research scale, we set the grid size to 1 km × 1 km in the regression model but neglected the fact that some landscape indices are sensitive to resolution. Furthermore, the threat factor in the InVEST model is subjective because of the difficulty in obtaining direct observation data; thus, there is room for improvement in the applied model.

The impact of urbanization on the environment is complex and variable. It is essential to focus on LUCC in long time series and take the effects of various factors into consideration, so as to make more targeted decisions of land management. This requires further studies on how to select landscape indices and set parameter values in order to reflect evolution features more precisely in the study area.

4.3. Proposals for Optimizing Urban Land Management

Through our results, it can be seen that the rapid expansion of construction land in the GBA from 1980 to 2020 has squeezed numerous ecological spaces. The main threat to regional ecological security has come from urbanization development and increasingly intense human activities, making the ecological environment foundation more fragile. Coastal lands managed and developed in a sustainable manner can provide national and regional benefits, notably in terms of achieving the Sustainable Development Goals (SDGs) [69]. Therefore, it is essential to divide functional areas rationally, exploit land resources moderately, and harmonize the relationship among cultivated land, construction land, and ecological land; this aspect deserves more attention from long-term urban development decision makers. With these issues in mind, according to the specific development features of the corresponding cities, a few suggestions for land management in the GBA are proposed:

• As the pilot cities of eco-environmental protection and policy reform, Guangdong and Shenzhen are required to integrate scattered land resources and improve land utilization efficiency. From an urban planning perspective, governments should divide areas by function and specify land use types (e.g., residential land, industrial and mining land, ecological protection land). Furthermore, it is urgent to establish red lines for construction land increments and strengthen the penalties for illegal construction. Meanwhile, governments should enhance cooperation with Hong Kong and Macao in ecological protection, promote the construction of public Transport-Oriented Development (TOD), and encourage green roof construction by learning from foreign experience. The ultimate goal is to achieve intensive and efficient development of mega-cities in the future.
• In the four cities along the Pearl River, i.e., Foshan, Dongguan, Zhongshan, and Zhuhai, reclamation activities of tidal flats and sea are common. Local governments should pay attention to the protection of coastal tidal flats and promote post-cultivation management. As for inter-provincial cooperation, it should rely on neighboring cities, strive to enhance cross-administrative cooperation in land planning, and strengthen wetland protection.
• Forestland resources are widely distributed in Zhaoqing, Jiangmen, and Huizhou, acting as ecological barrier for the GBA. These cities should develop a scientific and operational protection system for their natural resources, e.g., through strengthening soil and water conservation. Sustainable rural construction also deserves more attention. Governments should amend the extensive development of rural areas, strengthen the construction of public facilities, and encourage population clustering in city centers.

5. Conclusions

The evolution characteristics of land use and habitat quality in the GBA from 1980 to 2020 were explored in this study using the InVEST model. Furthermore, the impact of urbanization on habitat quality was quantitatively investigated using the GWR model. The results indicated that the GBA has experienced frequent and dramatic LUCC, dominated by the rapid expansion of construction land and the continuous reduction of cultivated land. These phenomena were more obvious in the southeastern cities which underwent rapid economic development, represented by Shenzhen and Dongguan. The overall landscape pattern has been moving on a diversity and balance way, influenced by the change of land use patterns, and the degree of landscape fragmentation increased. From 1980 to 2020, the habitat quality in the GBA showed a decreasing trend, with the most significant change occurring in the period of 2000 to 2010. The level of habitat quality in the outer areas was higher than that in the central regions. Human activities likely accelerate the degradation of habitat quality. Terrain conditions determined the overall distribution pattern of habitat quality. The effects of construction land, population, and road construction on habitat quality displayed significant spatially heterogeneity. We suggest that local governments should improve land use efficiency and engage in rational planning of urban functional zones in order to limit further deterioration of habitat quality.