Driving Factors and Risk Assessment of Rainstorm Waterlogging in Urban Agglomeration Areas: A Case Study of the Guangdong-Hong Kong-Macao Greater Bay Area, China

Abstract

Understanding the driving factors and assessing the risk of rainstorm waterlogging are crucial in the sustainable development of urban agglomerations. Few studies have focused on rainstorm waterlogging at the scale of urban agglomeration areas. We used the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) of China as a case study. Kernel density estimation (KDE) and spatial autocorrelation analysis were applied to study the spatial distribution characteristics of rainstorm waterlogging spots during 2013–2017. A geographical detector (GD) and geographically weighted regression (GWR) were used to discuss the driving mechanism of rainstorm waterlogging by considering eight driving factors: impervious surface ratio (ISR), mean shape index of impervious surface (Shape_MN), aggregation index of impervious surface (AI), fractional vegetation cover (FVC), elevation, slope, river density, and river distance. The risk of rainstorm waterlogging was assessed using GWR based on principal component analysis (PCA). The results show that the spatial distribution of rainstorm waterlogging in the GBA has the characteristics of multicenter clustering. Land cover characteristic factors are the most important factors influencing rainstorm waterlogging in the GBA and most of the cities within the GBA. The rainstorm waterlogging density increases when ISR, Shape_MN, and AI increase, while it decreases when FVC, elevation, slope, and river distance increase. There is no obvious change rule between rainstorm waterlogging and river density. All of the driving factors enhance the impacts on rainstorm waterlogging through their interactions. The relationships between rainstorm waterlogging and the driving factors have obvious spatial differences because of the differences in the dominant factors affecting rainstorm waterlogging in different spatial positions. Furthermore, the result of the risk assessment of rainstorm waterlogging indicates that the southwest area of Guangzhou and the central area of Shenzhen have the highest risks of rainstorm waterlogging in GBA. These results may provide references for rainstorm waterlogging mitigation through urban renewal planning in urban agglomeration areas.

1. Introduction

With the development of economic globalization and regional economic integration, the urbanization level has led to a new stage of urban agglomeration, which is a highly developed spatial form of integrated cities instead of a single city [1,2], and the natural surfaces of urban areas have undergone a significant and dramatic transformation, usually as a result of decreases in the natural vegetation cover and increases in impervious surfaces, which prevents rainwater from infiltrating into the soil [3,4,5]. Coupled with global climate change, subsequent waterlogging following rainstorms has become a serious problem in urban agglomeration areas worldwide [6]; waterlogging badly affects human life and property, and causes economic losses and environmental pollution [7,8]. In recent years, many influential urban agglomeration areas in the world have been seriously hit by rainstorm waterlogging and flood disasters. For instance, during late August 2011, Hurricane Irene, which was a large and destructive tropical cyclone, affected the New York Bay Area of the United States, resulting in widespread destruction and at least 49 deaths [9]. The 2017 California floods affected the San Francisco Bay Area and caused $73 million in damage [10,11]. During May 2014, cities in the Guangdong-Hong Kong-Macao Greater Bay Area (GBA) of China suffered from urban waterlogging and flood disasters due to heavy rain, badly affecting the lives and property of residents and bringing severe challenges to the sustainable development of the economy and society [12,13]. Therefore, rainstorm waterlogging disasters have become an “urban agglomeration disease” that must be mitigated by government and urban agglomeration planning agencies in the planning and construction stages of urban development.

Enhancing the resistance of an area against rainstorm waterlogging disasters has become an urgent mission for the development and construction of urban agglomeration areas. While many studies have focused on independent cities or on parts of cities and have applied these scales successfully, there is a lack of research on rainstorm waterlogging disasters at the scale of urban agglomeration areas. Many strategies have been proposed, such as drainage system improvements [14], low impact development [15], green infrastructure [16], and impervious surface optimization [17], to mitigate the losses caused by rainstorm waterlogging at the city scale. The effectiveness of such strategies depends on understanding the driving factors and influencing mechanisms behind rainstorm waterlogging [18,19,20]. Rainstorm waterlogging is a complex natural and social phenomenon that is the result of the dual effects of natural factors and human activity factors [21]. From the perspective of natural factors, rainstorm waterlogging is mainly affected by meteorological and topographical conditions. Among these factors, frequent extreme weather caused by global climate change results in high-intensity local rainfall concentrations [22,23]. Moreover, urban centers are often flat and lack the slopes required for drainage; thus, they have relatively poor drainage capacities [24]. From the perspective of human activity factors, poor drainage infrastructure and land use change have significant impacts on rainstorm waterlogging. Among these factors, current urban drainage systems generally have problems such as low design standards, unreasonable construction, and insufficient maintenance [25,26]. Hence, the drainage capacities of the current systems are insufficient to meet the drainage demand when extreme precipitation events occur. Moreover, increases in impermeable surfaces reduce drainage ability, resulting in increases in the frequency and intensity of rainstorm waterlogging events [27,28]. Furthermore, woodland, grassland, and other green areas have significant effects on intercepting rainfall, reducing runoff and delaying confluence, and can effectively reduce the risk of waterlogging [29]. In addition to the abovementioned factors, changes in urban water area have a regulating effect on waterlogging, thereby affecting the occurrence of rainstorm waterlogging events [30].

Various studies have been conducted on the driving factors of rainstorm waterlogging in typical cities that have been impacted by rainstorm waterlogging and flood disasters for a long time. Tehrany et al. [31] indicated that elevation and slope have significant impacts on rainstorm waterlogging. Gaitan et al. [32] confirmed that, in addition to the significant influences of elevation and slope on rainstorm waterlogging, land use/land cover is another key factor that drives rainstorm waterlogging. Land use/land cover has an even greater impact than elevation and slope [33]. With rapid urbanization, the impacts of land use/land cover changes, and the rapid expansion of impervious surfaces are particularly significant. Land use/land cover changes have gradually become the main reasons for increasingly serious rainstorm waterlogging disasters [34]. Su et al. [35] explored the cause of rainstorm waterlogging in Kunming and found that the increase in impervious surfaces, such as construction land, will aggravate the risk of rainstorm waterlogging due to the impermeability of the construction land and the rapid speed of rainstorm runoff. Furthermore, there are also some studies focusing on related work, i.e., the impact of the spatial pattern of impervious surfaces on the occurrence of rainstorm waterlogging. These studies analyze how the area, proportion, and spatial distribution pattern of impervious surfaces affect surface runoff and rainstorm waterlogging disasters [36,37]. Zhang et al. [38] took Guangzhou as a case study and used redundancy analysis to investigate the influence of the spatial pattern of impervious surfaces on rainstorm waterlogging. The results revealed that the composition of the impervious surfaces contributed more to rainstorm waterlogging than did the configuration of the total impervious surface. Yu et al. [39] also explored the spatiotemporal relationship between impervious surface expansion and urban rainstorm waterlogging by using local spatial autocorrelation (LSA) and geographically weighted regression (GWR) in downtown Guangzhou, and found that the impacts of impervious surfaces on urban waterlogging were quite different among different subregions of the research area. Therefore, there is spatial non-stationarity in the relationship between driving factors and rainstorm waterlogging, especially in regions with complicated geographical environments. Wu et al. [40] indicated the different impacts of different regional landscape patterns on rainstorm waterlogging in Shenzhen by using the GWR model, providing useful information for urban planning and mitigating waterlogging disasters.

Generally, previous studies have described the influences of different factors on rainstorm waterlogging and the spatial non-stationary nature of these influences at the city scale, but this may not be sufficient considering the complex environment of urban agglomerations. The occurrence of rainstorm waterlogging is not limited to any independent city, but is always concentrated in a large-scale region, such as urban agglomeration and metropolitan areas, which have much more complicated environmental features. A few studies have explored urban expansion and the evolution of ecological quality at the urban agglomeration scale and found that there are large differences among subregions within the large-scale region [5,41,42,43,44]. Therefore, it will be beneficial to study the driving factors and influencing mechanisms of rainstorm waterlogging in urban agglomeration areas, and the results will be helpful for predicting rainstorm waterlogging areas in which an effective way to reduce or even prevent waterlogging-related losses can be provided.

The GBA is the largest morphologically contiguous urban agglomeration in the world; the area has complex spatial patterns and topographical features, and it is constructed to enable Guangdong, Hong Kong, and Macao to achieve complementary advantages, coordinated development, and mutual benefits [41,45,46,47]. Since the formal implementation of China’s reform and opening up in 1978, the GBA has experienced remarkable urbanization and industrialization, resulting in a dramatic change in land use/land cover from vegetative cover to urban areas [48,49]. It is reported that rainstorm waterlogging events abound in the GBA, and that these events badly influence the lives and property of residents and introduce severe challenges to the sustainable development of the economy and society [12,50,51]. Therefore, the GBA is a typical case that can be studied in order to reach a better understanding of the driving factors and influencing mechanisms of rainstorm waterlogging at the scale of the urban agglomeration area. Thus, in this study, kernel density estimation (KDE) and spatial autocorrelation were utilized to explore the spatial characteristics of urban rainstorm waterlogging in the GBA. Then, the degree of influence of each driving factor on rainstorm waterlogging and the interactions between the factors were discussed by using the geographical detector (GD). Finally, the geographically weighted regression (GWR) model was used to reflect the non-stationarity of the driving factors impacting urban rainstorm waterlogging and to assess the risk of rainstorm waterlogging in the GBA based on principal component analysis (PCA). This study aims to explain the following issues:

(1)

What are the main driving factors of rainstorm waterlogging in the GBA? How do the various driving factors affect the occurrence of rainstorm waterlogging?

(2)

Are there spatial differences in the impacts of various driving factors on rainstorm waterlogging in the GBA?

(3)

Where are the areas with high risks of rainstorm waterlogging in the GBA under the combined effect of various driving factors?

(4)

How can the risk of rainstorm waterlogging in the GBA be effectively alleviated?

2. Study Area and Datasets

2.1. Study Area

The GBA is the most dynamic and urbanized urban agglomeration in China, and includes 11 cities (Dongguan, Foshan, Guangzhou, Huizhou, Jiangmen, Shenzhen, Zhaoqing, Zhongshan, Zhuhai, Hong Kong, and Macao) (Figure 1) [47]. It is located between 21°30′~24°30′ N and 111°15′~115°30′ E and covers approximately 56,000 km². The terrain of the GBA decreases in elevation from north to south, with hills in the north, east, and west and flood plains in the central region and in the south. Moreover, the GBA is located in the lower reaches of the East River, North River, and West River in the Pearl River system, with complicated drainages and densely covered waterways. The GBA is in an area with a tropical and subtropical monsoon climate and plenty of rainfall. The annual average precipitation in the GBA is approximately 1500 mm, and heavy precipitation occurs from June to October under the effects of typhoons and tropical pressures [52]. Since the reform and development of China, the urbanization of the GBA has accelerated. Coupled with the occurrence of regional short-term extreme rainfall events, most of the cities in the GBA have frequently suffered from rainstorm waterlogging, causing substantial losses of human life and property under the backdrops of climate change and rapid urbanization [53,54].

2.2. Datasets and Pre-Processing

The acquisition of multi-source data is key to better understanding and evaluating rainstorm waterlogging risk in the GBA. In this study, the raw data sources utilized contain remote sensing images, digital elevation maps, text information of waterlogging events, etc. All raw data sources are shown in

Table 1. Metadata information.

Data	Data Type	Data Year	Data Source
Waterlogging information	Text file	2013–2017	City water authority, Mainstream media
Landsat remote sensing image	Satellite imagery	2017	USGS (https://www.usgs.gov/) (accessed on 25 January 2021)
High-resolution remote sensing image	Satellite imagery	2017	Google Earth satellite imagery
Topographic data	Raster	2009	Geospatial Data Cloud (http://www.gscloud.cn/) (accessed on 25 January 2021)
River data	Vector	2017	Open Street Map (https://www.openstreetmap.org/) (accessed on 25 January 2021)
1:4 million administrative divisions	Vector	2005	National Geomatics Center of China (http://www.ngcc.cn/ngcc/) (accessed on 25 January 2021)

.

2.2.1. Rainstorm Waterlogging Inventory

In this study, we obtained rainstorm waterlogging records with accurate time and address information from the news media and water authorities in cities within the GBA. We used ArcGIS and Google Earth to visually identify and map the spatial locations of these rainstorm waterlogging areas. In the process of mapping these rainstorm waterlogging spots, if the address of the waterlogging spot was described as a specific road intersection, community entrance, concrete building, etc., we represented the waterlogging spot with a dot. If the address of the waterlogging spot was a road section, we regarded the geometric center of the road as the spatial location of the waterlogging spot. Approximately 954 rainstorm waterlogging spots from 2013–2017 were found (Figure 2).

2.2.2. Driving Factors of Rainstorm Waterlogging

Considerable research has indicated that factors driving the frequent occurrence of rainstorm waterlogging events can be summarized as follows: (1) land use/land cover change caused by urbanization [55]; (2) spatial variation of terrain surface relief [56]; (3) degradation of the drainage function of the river network [57]; and (4) insufficient drainage capacity of the underground pipeline network [58]. Eight driving factors were selected in this study after a review of previous studies; these include the impervious surface ratio (ISR), mean shape index of the impervious surface (Shape_MN), aggregation index of impervious surface (AI), fractional vegetation cover (FVC), elevation, slope, river density, and river distance. These driving factors can be categorized into three classes: land cover characteristics (composition and spatial configuration), topography, and water network.

Land Cover Characteristics

During the process of rapid urbanization, land use/land cover types have changed considerably, and the original natural or semi-natural land use/land cover types have been replaced by hard surfaces. Such a transformation has a destructive effect on the natural hydrologic cycle by decreasing the surface permeability, which may increase the risk of rainstorm waterlogging [54]. Impervious surfaces are typical hard surfaces through which rainwater infiltration is much lower than that through pervious surfaces [59]. Several studies have shown that both the composition and configuration of impervious surfaces reduce the infiltration of water into the surface, and thereby increase the associated risk of rainstorm waterlogging [60,61]. Therefore, the spatial pattern of impervious surfaces was designated as a driving factor. In this paper, the modified linear spectral mixture analysis (MLSMA) method was used to extract the impervious surface data of the GBA by utilizing Landsat data from 2017. Then, the impervious surface mapping accuracy was assessed after comparison with high-resolution data acquired from Google Earth. The accuracy of the classification result was 88.3%, which indicated that the accuracy of the impervious surface extraction by the MLSMA was sufficient for further analysis. The landscape pattern index can be used to quantitatively study the composition and configuration of the landscape. Hence, many studies have utilized landscape pattern indexes to extensively explore the spatial patterns of impervious surfaces [62,63]. Considering the actual connotation of the landscape pattern and actual impact of the impervious surface spatial pattern on rainstorm waterlogging, ISR, Shape_MN, and AI were selected as the driving factors of rainstorm waterlogging based on the principle of mutual independence of the factors representing landscape patterns [64]. The implication of these indices is listed in

Table 2. Overview of indices used in the current study.

Name	Range	Implication
ISR	1 ≥ ISR ≥ 0	ISR is the ratio of the impervious surface. Bigger values mean a higher impervious surface area.
Shape_MN	Shape_MN ≥ 1	Shape_MN represents the mean shape value of the landscape. The larger Shape_MN is, the more irregular the landscape distribution.
AI	100 ≥ AI > 0	AI describes the aggregation degree of landscape, and helps characterize the landscape in terms of stability. Higher values mean greater clustering.

. Among them, ISR was obtained based on the class area (CA), which was calculated by Fragstats [65]. Both Shape_MN and AI were calculated by using Fragstats as well.

Green space is another typical land use/land cover type that can represent surface permeability, and has a negative impact on rainstorm waterlogging. However, due to rapid development, green space has been continuously reduced, which in turn affects the infiltration of surface rainwater, resulting in an increased risk of rainstorm waterlogging. In this paper, the FVC was used to reflect green space coverage in the GBA. The FVC is defined as the percentage of vegetation occupying a unit area, and is calculated based on the normalized difference vegetation index (NDVI) derived from remote sensing data [66]. The FVC is calculated as follows:

$$\mathrm{FVC}=\frac{{\mathrm{NDVI}}_i-{\mathrm{NDVI}}_{min}}{{\mathrm{NDVI}}_{max}-{\mathrm{NDVI}}_{min}}$$

(1)

where ${\mathrm{NDVI}}_i$ is the NDVI value of the $i$th pixel, ${\mathrm{NDVI}}_{max}$ and ${\mathrm{NDVI}}_{min}$ are the maximum and minimum values of NDVI, respectively, and the range of the FVC is [0, 1].

Topography

Elevation is considered to be one of the most influential factors affecting rainstorm waterlogging [67]. In general, rainstorm waterlogging events are more likely to occur in regions with relatively lower elevations, because these regions easily obtain surface runoff [31]. In addition, the slope is also an important topographical feature that impacts the occurrence of rainstorm waterlogging [68]. A region with a steeper slope has a lower possibility of being inundated because rainwater can easily flow downslope [69]. Therefore, in this study, we selected elevation and slope to represent the topographic characteristics. The topographic factors in this study were derived from ASTER GDEM, and have resolutions of 30 m.

Water Network

The water network structure is related to the ability of an area to store rainwater storage and to the ecological carrying capacity of an area. Many rivers, canals, and lakes have been converted into built-up areas due to urban construction and human activities; this conversion destroys the water network structure and decreases the water storage capacity. Therefore, the water network distribution is closely related to rainstorm waterlogging. River density and river distance are common structural indexes used to characterize the water storage capacity of rivers [57]. In this study, water network factors were obtained based on the Open Street Map Dataset. The river density can be measured as follows:

$$D_r=\frac{L}{F}$$

(2)

where $D_r$ represents the river density in km/km², $L$ is the total river length in km, and $F$ is the total area in km². We calculated the river distance, which is the distance to the river in the GBA, by using the Euclidean distance method.

3. Methods

3.1. Integrated Frameworks

In this paper, the complex relationships between rainstorm waterlogging and the driving factors of waterlogging, as well as the spatial variability of these relationships, were explored. The integrated framework of this research is shown in Figure 3. First, KDE and spatial autocorrelation analysis were utilized to explore the spatial distribution characteristics of rainstorm waterlogging in the GBA. Second, eight driving factors (i.e., ISR, Shape_MN, AI, FVC, elevation, slope, river density, and river distance), which can be divided into three classes (i.e., land cover characteristics, topography, and water network), were selected in this study by reviewing previous studies to investigate the driving mechanisms of rainstorm waterlogging. Then, GD was carried out to reveal the relative contributions of and the interactions between the driving factors. Finally, GWR was carried out to reflect the spatial non-stationarity of the relationships between rainstorm waterlogging and the driving factors, and to assess the spatial distribution of rainstorm waterlogging risk in the GBA. More details are given in the following sections.

3.2. Spatial Pattern Analysis Method

Rainstorm waterlogging is influenced by the urban environment and by human activities; hence, it has the notable characteristic of aggregated distribution in local regions. To provide basic information for preventing rainstorm waterlogging and encouraging urban renewal, spatial analysis methods were utilized in this study to explore the spatial pattern of waterlogging spots. Because the waterlogging spots were handled as points, it was impossible to visually express the waterlogging situation around the waterlogging spots. Therefore, KDE [70] was used to obtain the rainstorm waterlogging degree in the GBA, which can better reflect the impact of a rainstorm waterlogging event in the vicinity.

The KDE method is a density distribution model with the spatial feature of the distance attenuation effect, which can be used to compute the density of randomly distributed points in a neighborhood around each output raster pixel. It is functioned by setting the search scope (i.e., window), which is a circular area with a radius of bandwidth around the raster pixel in the study area. Based on the principle of anti-distance weight, the central raster pixel of the window gives the weight of each raster pixel outwards. Then, the kernel density value of the central raster pixel is obtained by summing the product of all of the kernel density values and weights in the window. Considering the applicability of the study, the Rosenblatt–Parzen estimation method, which is the most common KDE, was applied in this paper; its formula is as follows:

$$f_n\left(X\right)=\frac{1}{nh}{{\displaystyle \sum }}_{i=1}^{n}K\left(\frac{X-X_i}{h}\right)\left(i=1,2,\dots \dots n\right)$$

(3)

Let $f_n\left(X\right)$ be the estimated kernel density of the central raster pixel $X$, where $n$ is the number of waterlogging spots in the window, $h$ is the bandwidth, $K$ is the kernel function, and $\left(X-X_i\right)$ is the distance from the central raster pixel $X$ of the window to the waterlogging spot raster pixel $X_i$. In KDE, the bandwidth is an important parameter. If the bandwidth is too large, the point density surface will become too smooth to cover the density structure of the point pattern, while a bandwidth that is too small will cause the point density distribution to change very abruptly [71]. Therefore, the optimal bandwidth is determined by repeatedly setting the bandwidth and comparing the smoothness of the point density surfaces. The kernel function is a quadratic kernel function, and the formula is as follows:

$$K\left(X\right)=\{\begin{array}{c}\frac{3}{4}\left(1-X^2\right), X\in \left[0, 1\right]\\ 0 , otherwise\end{array}$$

(4)

The distribution of waterlogging spots in the GBA is concentrated as shown in Figure 2, which implies that the occurrence of waterlogging events has a certain dependence on space. The dependence of waterlogging spots on space and its spreading characteristics can be understood by visualizing the spatial agglomeration degree of rainstorm waterlogging. Therefore, we meshed the study area into grids with a cell size of 1 km in this study. After assigning the KDE result to the grids in the study area, we used the spatial autocorrelation analysis method to measure the spatial dependence of waterlogging spots, which includes global spatial autocorrelation (GSA) and local spatial autocorrelation (LSA). GSA is utilized to generalize the overall degree of spatial dependence for waterlogging spots in the whole study area. Global Moran’s I is an important index to measure GSA. The calculation method is as follows:

$$I=\frac{{\sum }_{i=1}^n{\sum }_{j=1}^n{w}_{ij}\left(y_i-\overline{y}\right)\left(y_j-\overline{y}\right)}{S^2\left({\sum }_i{\sum }_j{w}_{ij}\right)}$$

(5)

In this formula, $n$ is the number of grid cells indexed by $i$ and $j$ across the global, and $y_i$ and $y_j$ are the kernel density of rainstorm waterlogging spots in grid cells $i$ and $j$, respectively. $\overline{y}$ is the mean value of the kernel density of rainstorm waterlogging. $w_{ij}$ is the element of the weight matrix $w$, indicating the spatial weights of the grid cells $i$ and $j$, where $w_{ij }$= 1 if grid cells $i$ and $j$ are adjacent and $w_{ij}$= 0 otherwise. The adjacent relation is of the Queen type, that is, grid cells are adjacent if they share a common border or vertex [72]. Meanwhile, $S^2$ is the discrete variance of $y_i$, which can be calculated using Equation (6):

$$S^2=\frac{1}{n}{{\displaystyle \sum }}_{i=1}^n{\left(y_i-\overline{y}\right)}^2$$

(6)

The Z-test value was used to measure the significance of Moran’s I using Equation (7):

$$Z=\frac{I-E\left(I\right)}{\sqrt{\mathrm{Var}\left(I\right)}}$$

(7)

where $E\left(I\right)$ represents the mathematical expectation of Moran’s I under the hypothesis of spatial stochastic distribution and $\mathrm{Var}\left(I\right)$ is the variation coefficient.

LSA is utilized to estimate the variation of the autocorrelation degree for rainstorm waterlogging spots across a geographic space. Local Moran’s I is used as a local indicator of spatial association (LISA) to reveal the degree of LSA for waterlogging spots by evaluating the similarity of rainstorm waterlogging density for each grid cell with that of its surroundings [73,74]. The calculation method is as follows:

$$I_i=\frac{\left(y_i-\overline{y}\right){\sum }_{j=1}^n{w}_{ij}\left(y_j-\overline{y}\right)}{S^2}$$

(8)

Generally, four spatial cluster/outlier types of rainstorm waterlogging spots are classified using the local Moran’s I, including the high-high cluster, low-low cluster, high-low outlier, and low-high outlier. Among them, the high-high cluster indicates that the location has high values of the waterlogging spot kernel density and a high level of similarity with its surroundings; the low-low cluster indicates that the location has low values of the waterlogging spot kernel density and a high level of similarity with its surroundings; the high-low cluster indicates that the location has high values of the waterlogging spot kernel density and a low level of similarity with its surroundings; and the low-high cluster indicates that the location has low values of the phenomenon and a low level of similarity with its surroundings. Therefore, high-high and low-low clusters reflect positive spatial correlation, while high-low and low-high outliers reflect negative spatial correlation. The spatial cluster/outlier types of different locations in the study area can be represented by a LISA cluster map [75].

3.3. Mechanism Analysis Method

The causes of rainstorm waterlogging events are often multiple and intricate. Studying the specific causes of rainstorm waterlogging events and the interactions between them is helpful for providing reasonable suggestions for the prevention of rainstorm waterlogging and the optimization of urban renewal in the future. In addition, there are dependencies among the driving factors of rainstorm waterlogging. Hence, we utilized a GD that is not affected by the collinearity of factors to quantitatively assess the impacts of driving factors on rainstorm waterlogging [76].

GD is a new statistical method used to detect spatially stratified heterogeneity and reveal the driving factors behind the heterogeneity; it includes a risk factor detector, a risk area detector, an ecological detector, and an interaction detector [77,78]. In this study, the risk factor detector, risk area detector, and interaction detector of the GD were selected. The risk factor detector can be used to identify the dominant factors influencing rainstorm waterlogging in the study area, and the formula is as follows:

$$q=1-\frac{{\sum }_{h=1}^L{N}_h{\sigma }_h^2}{N{\sigma }^2}$$

(9)

In this formula, $q$ represents the contribution of the driving factor (i.e., ISR, Shape_MN, AI, FVC, elevation, slope, river density, and river distance) to rainstorm waterlogging with the value range [0, 1], and the greater the value of $q$, the greater the contribution to rainstorm waterlogging. Due to the requirement of GD that the input variable should be sub-regions (statistical units), the driving factor should be stratified into $h=1, 2\dots L$ strata. In this paper, the continuous value of each driving factor were discretized into nine strata by using the natural discontinuity (ND) method referencing the experience of optimal discretization for GD, which was proposed by Wang and Xu [76]. Moreover, $N$ is the total number of grid cells in study area, while $N_h$ denotes the number of grid cells in stratum $h$. ${\sigma }^2$ and ${\sigma }_h^2$ denote the global variance and the divisional variations of the dependent variable (i.e., the kernel density of rainstorm waterlogging), respectively.

A risk area detector can be used to determine whether different strata of factors have variable impacts on the spatial distribution of rainstorm waterlogging by using a t-test [79]. The formula is as follows:

$$t_{{\overline{y}}_{h=i}-{\overline{y}}_{h=j}}=\frac{{\overline{Y}}_{h=i}-{\overline{Y}}_{h=j}}{{\left[\frac{\mathrm{Var}\left({\overline{Y}}_{h=i}\right)}{n_{h=i}}+\frac{\mathrm{Var}\left({\overline{Y}}_{h=j}\right)}{n_{h=j}}\right]}^{\frac{1}{2}}} \left(i=j=1,2,\dots \dots L\right)$$

(10)

where ${\overline{Y}}_h$ represents the mean value of the rainstorm waterlogging kernel density in stratum $h$, $n_h$ is the number of samples in stratum $h$, and $\mathrm{Var}$ represents the variance. To analyze how the driving factors affect rainstorm waterlogging, the mean values of the waterlogging degree in different strata (i.e., ${\overline{Y}}_h$) of each driving factor were obtained by using the risk area detector in this paper.

To evaluate the impact of the interaction between different driving factors on rainstorm waterlogging, an interaction detector was used in this study. The processing can be carried out as follows: first, the $q$ value of two driving factors $X_1$ and $X_2$ to $Y$ (i.e., $q\left(X_1\right)$ and $q\left(X_2\right)$) are calculated, respectively, then, the $q$ value of the interactions between the driving factors (i.e., $q\left(X_1{\displaystyle \cap }{X}_2\right)$) is operated by merging the layers of the two driving factors $X_1$ and $X_2$ via the previous formula. The interaction relationships between the two factors are as follows: if $q\left(X_1{\displaystyle \cap }{X}_2\right)<\mathrm{Min}\left(q\left(X_1\right),q\left(X_2\right)\right)$, it is shown as nonlinear weakening, in which case the conclusion is that the contributions of joint $X_1$ and $X_2$ non-linearly weaken the contribution of each another; if $\mathrm{Min}\left(q\left(X_1\right),q\left(X_2\right)\right)<q\left(X_1{\displaystyle \cap }{X}_2\right)<\mathrm{Max}\left(q\left(X_1\right),q\left(X_2\right)\right)$, it is shown as unique weakening, which indicates that the $X_1$ and $X_2$ joint contribution non-linearly weakens the single contribution; if $q\left(X_1{\displaystyle \cap }{X}_2\right)>\mathrm{Max}\left(q\left(X_1\right),q\left(X_2\right)\right)$, it is shown as binary enhancing, in which case the conclusion is that the contributions of joint $X_1$ and $X_2$ enhance the contribution of each another; if $q\left(X_1{\displaystyle \cap }{X}_2\right)= q\left(X_1\right)+q\left(X_2\right)$, it is shown as independent, which indicates that the $X_1$ and $X_2$ joint contribution doesn’t affect the single contribution; and if $q\left(X_1{\displaystyle \cap }{X}_2\right)> q\left(X_1\right)+q\left(X_2\right)$, it is shown as nonlinear enhancing, in which case the conclusion is that the $X_1$ and $X_2$ joint contributions non-linearly enhance the contribution of each another.

In this study, eight driving factors were selected from the three classes (i.e., land cover characteristics, topography, and water network) as the independent variables; the differences in the contributions of each of these factors to rainstorm waterlogging occurrence in urban agglomeration areas and the interactions between the factors were detected.

3.4. Regression Analysis Method

Due to the differences in the development factors of cities in the GBA, the spatial distribution of both rainstorm waterlogging spots and driving factors are characterized by non-stationary features; that is, the correlative degree between the driving factors and the occurrence of rainstorm waterlogging varies with the change in spatial location. Therefore, in this paper, based on the analysis of the driving factors of rainstorm waterlogging obtained by using the GD, we also carried out GWR analysis and explored the spatial distribution of rainstorm waterlogging risk in the GBA. GWR is a modelling technology that can effectively deal with geospatial differences, and is a spatial extension of the ordinary least squares regression (OLS) model [80]. GWR returns the spatial position information to the parameters to estimate the local relationships, thereby obtaining a more objective result [81]. Hence, the GWR model is an effective tool for studying the driving factors of rainstorm waterlogging in different geographical locations. The model can be expressed as follows:

$$y_i={\beta }_{i0}+{\displaystyle \sum }_{k=1}^m{\beta }_{ik}{x}_{ik}+{\epsilon }_i\left(i=1,2,3\dots \dots n\right)$$

(11)

where $y_i$ represents the dependent variable, which is the kernel density of rainstorm waterlogging in grid cell $i$ in this paper, ${\beta }_{i0}$ is the intercept of the grid cell $i$, $m$ represents the total number of explanatory variables, which are driving factors of rainstorm waterlogging, ${\beta }_{ik}$ represents the regression coefficient of the $k$th explanatory variable in the grid cell $i$, $x_{ik}$ represents the $k$th explanatory variable of the grid cell $i$, and ${\epsilon }_i$ represents the random error of the grid cell $i$. The spatial correlations between the model residuals were judged by the Z-test value, and the model could be trusted when the residuals were randomly distributed. Spatial weighting is the key to GWR, and its functions are mainly divided into fixed bandwidth and adaptive weight functions. The adaptive weight function is adopted in this study. The formula is as follows:

$$w_{ij}=\{\begin{array}{c}1, j\in k\\ 0, j\notin k\end{array}$$

(12)

where $k$ is the number of elements abutting element $i$. In this study, the optimal number of adjacent elements was set according to the minimum AICC principle proposed by Fotheringham [80].

However, when GWR is used for evaluation, more input variables lead to high computational complexity, while fewer input variables lead to low accuracy of the model. Moreover, there are often certain correlations between the driving factors of rainstorm waterlogging. Li [82] proposed the GWR model based on PCA, which not only eliminates the correlations between the geographical elements, but also improves the explanatory ability of the model. PCA is a common method to solve the multicollinearity problem in multiple linear regression. It is a method of integrating multiple related variables into a few unrelated variables based on the principle of as little information loss as possible. The PCA-GWR method combines the advantages of PCA and GWR. Before carrying out GWR to do a risk assessment of rainstorm waterlogging, PCA was used to generate the unrelated principal components by reducing the dimensionality of the variables, which were driving factors.

Therefore, this paper first used the GWR model to study the spatiotemporal differentiation of the impacts of driving factors on rainstorm waterlogging. Then, the PCA-GWR method was used to assess the risk of rainstorm waterlogging in the GBA by carrying out GWR.

4. Result

4.1. Spatial Distribution Characteristics of Rainstorm Waterlogging Spots in the GBA

The total number of rainstorm waterlogging spots was 954 in the GBA during the period from 2013–2017. The number of rainstorm waterlogging spots in Shenzhen was the highest, which had 280 rainstorm waterlogging spots, followed by Guangzhou, which had 171 rainstorm waterlogging spots. Jiangmen, Foshan, and Zhuhai also had a large number of rainstorm waterlogging spots, at 108, 97, and 97 spots, respectively. The numbers of rainstorm waterlogging spots in Dongguan and Huizhou were relatively small, at 67 and 62 spots, respectively. There were only 39, 22, and 12 rainstorm waterlogging spots in Zhaoqing, Zhongshan, and Hong Kong, respectively. Macao had no rainstorm waterlogging spots, according to the statistics used in this study.

The kernel density distribution of rainstorm waterlogging is shown in Figure 4a. It shows that the rainstorm waterlogging density in the GBA presents a multicore and multilevel spatial aggregation structure. The highest rainstorm waterlogging density is mainly concentrated in the southwest portion of the downtown area of Guangzhou and in the middle of Shenzhen. Two core areas with relatively higher densities of rainstorm waterlogging are formed in the connecting area of Guangzhou and Foshan, the whole city of Shenzhen, and the adjacent area of Shenzhen and Hong Kong. Meanwhile, three secondary core areas with high density values are formed in northern Dongguan, northeastern Jiangmen, and eastern Zhuhai. The northwest region of Foshan, the southwest area of Huizhou, the southeast area of Zhaoqing, and the west region of Zhuhai form transitional areas with medium densities of rainstorm waterlogging. Other regions are peripheral areas with low densities of rainstorm waterlogging.

The Global Moran’s I and Z-test value were 0.995 and 809.45, respectively, which reveal that the occurrence of rainstorm waterlogging in the GBA was distributed in a clustering spatial form. In addition, after calculating the local Moran’s I, we found that there were two cluster types of the waterlogging spots pattern in GBA, which were the high-high cluster and low-low cluster. To analyze the spatial aggregation of high waterlogging degree area in the GBA, the high-high cluster result was exported and mapped as shown in Figure 4b. It shows that there are three high-high cluster areas with the distribution characteristics of large-area gathering and axial extension. The first area is located in the connecting area of Guangzhou and Foshan, the second is located in the south of Zhongshan and Zhuhai, and the third is located in southern Dongguan, southwest Huizhou, north Hong Kong, and all of Shenzhen. There are also several relatively small circular high-high cluster areas distributed in various cities in the GBA. The above results indicate that the distribution of rainstorm waterlogging hot spots in the GBA is concentrated.

4.2. Factors Driving Rainstorm Waterlogging in the GBA

In this study, the GD was used to explore the effect intensities and the interactive effects of the driving factors on rainstorm waterlogging in both the whole GBA and the cities in the GBA. It should be noted that since there was no rainstorm waterlogging event record in Macao during the study period, it is impossible to explore the mechanisms impacting rainstorm waterlogging in Macao.

4.2.1. Dominant Driving Factor Analysis

The influence of each driving factor on rainstorm waterlogging in the whole GBA and the cities within the GBA was calculated by the risk factor detection. The results are shown in

Table 3. The impact of driving factors on rainstorm waterlogging in the whole GBA and local cities.

	GBA	Dongguan	Foshan	Guangzhou	Hong Kong	Huizhou
ISR	0.1141	0.0328	0.1601	0.1769	0.0313	0.1780
Shape_MN	0.0472	0.0120	0.0514	0.0758	0.0225	0.1028
AI	0.0013	0.0214	0.1046	0.0253	0.0093	0.1525
FVC	0.0103	0.0242	0.1085	0.1381	0.0861	0.1901
Elevation	0.0022	0.0352	0.0666	0.0332	0.0536	0.0198
Slope	0.0029	0.0712	0.0859	0.0573	0.1158	0.0147
River Density	0.0064	0.0025	0.0211	0.0203	0.0606	0.0038
River Distance	0.0234	0.0035	0.0536	0.0016	0.0595	0.0236
Variation
	Jiangmen	Shenzhen	Zhaoqing	Zhongshan	Zhuhai
ISR	0.1774	0.0823	0.0735	0.0050	0.0224
Shape_MN	0.0724	0.0560	0.0402	0.0053	0.0082
AI	0.0326	0.0042	0.0024	0.0067	0.0057
FVC	0.3025	0.2021	0.1289	0.0613	0.0251
Elevation	0.0200	0.0617	0.0634	0.0041	0.0273
Slope	0.0207	0.1066	0.0434	0.0082	0.0459
River Density	0.0024	0.0043	0.0099	0.0016	0.0246
River Distance	0.0054	0.0595	0.0331	0.0073	0.0721
Variation

, with mini-graphs representing the changing trend of all of the data in the regions. As shown in the overall mini-graphs, the ranking of the effect of driving factors on the rainstorm waterlogging in the whole GBA and all of the cities within the GBA are totally different. For the whole GBA, the land cover characteristic factors, except AI, are determined to be the primary factors affecting the occurrence of rainstorm waterlogging. The secondary factors are the water network factors, and the explanatory power of river distance is much greater than that of river density. The explanatory powers of the topographic factors are relatively small, and the explanatory power of slope is greater than that of elevation. Land cover characteristic factors are the primary factors affecting the occurrence of rainstorm waterlogging in Foshan, Guangzhou, Huizhou, Jiangmen, and Zhaoqing, with ISR and FVC having the largest explanatory powers. The secondary factors are topographic factors. The explanatory power of slope in Foshan and Guangzhou is significantly higher than that of elevation, while the explanatory powers of elevation and slope in Huizhou and Jiangmen are not obvious. The explanatory power of elevation in Zhaoqing is significantly higher than that of slope. Water network factors have the least influence on the spatial distribution of rainstorm waterlogging in the above five cities. Both land cover characteristic factors and topographic factors in Dongguan and Shenzhen, located on the east bank of the Pearl River Estuary, are the primary factors affecting the spatial distribution of rainstorm waterlogging, while the water network factors have limited explanatory power. However, for Zhongshan and Zhuhai, which are located on the west bank of the Pearl River, the explanatory power of river distance is relatively large, and is even larger than that of the land cover characteristic factors. The analysis of the relative importance of the driving factors affecting rainstorm waterlogging in Hong Kong shows no obvious pattern. ISR, Shape_MN, and AI all have little explanatory power for the spatial distribution of rainstorm waterlogging, which is very different from with the results of other cities.

To further detect how each driving factor affects the spatial distribution of rainstorm waterlogging, the risk area detector of the GD was used. The results are shown in Figure 5. In this figure, a different color of the bar represents a different stratum, in which the larger the stratum value is, the greater the value of the driving factor. The vertical axis is the mean value of the rainstorm waterlogging degree in different stratum, while the horizontal axis represents different regions. Meanwhile, it is necessary to explain that the effect in some strata of driving factors may be interfered with by other driving factors, which may increase or decrease the rainstorm waterlogging degree in these strata, whereas the overall trend of all strata will not be influenced. Hence, the risk area detector results should be used to assist in the analysis of the driving mechanisms of the factors. The results show that there are significant differences in the degrees of rainstorm waterlogging between the different strata of each driving factor, which have some change rules. For the land cover characteristic factors, the degree of rainstorm waterlogging in the whole GBA and in most of the cities in the GBA increases with increasing ISR, Shape_MN, and AI. The lower the FVC is, the higher the possibility of rainstorm waterlogging events. For topographic factors, the areas with low elevation and low slope have high degrees of rainstorm waterlogging due to heavy rain. For water network factors, there is no obvious change rule between the degree of rainstorm waterlogging and river density. The closer the river distance is, the higher the possibility of rainstorm waterlogging occurrence.

4.2.2. Interactive Analysis of Driving Factors

The interaction detector was used in this study to examine the interactions between the driving factors of rainstorm waterlogging. The results are shown in Figure 6. The results reveal that the combination of different factors exerts a greater influence than does a single factor in the whole GBA and in each city within the GBA, and there are two main types of interactions found in this study: binary enhancing and nonlinear enhancing interactions. For the interaction types, with the exception of the whole GBA, Foshan, Guangzhou, and Zhaoqing, binary enhancement mainly occurs within each category of driving factors, such as the interaction between elevation and slope and between river density and river distance. In the whole GBA and in Guangzhou, binary enhancement also appears between land cover characteristic factors and water network factors, while in Foshan and Zhaoqing, binary enhancement appears among all three driving factor classes.

The intensities of the interactions between the driving factors of rainstorm waterlogging in the whole GBA and in the cities in the GBA are divided into three classes according to the principle of equal spacing, representing weak, moderate, and strong impacts of the driving factors, respectively. For the GBA, the interactions between ISR and other driving factors are the strongest, the Shape_MN of impervious surfaces has a strong interaction with the other driving factors, and the interactions between other driving factors are relatively weak, which reveal that the area and spatial distribution of impervious surfaces play decisive roles in influencing rainstorm waterlogging in the GBA. The intensities of the interactions between driving factors have obvious differences in the cities within the GBA. For Dongguan, Foshan, Guangzhou, Shenzhen, and Zhongshan, strong interactions mainly occur between land cover characteristic factors and topographic factors. Among these factors, the interaction between ISR and slope is the strongest in Dongguan and Foshan, while the interaction between FVC and slope is the strongest in Guangzhou and Zhongshan; in Shenzhen, the interaction between FVC and slope is the strongest. In Huizhou, Jiangmen, and Zhaoqing, strong interactions mainly occur between ISR, FVC, and other driving factors, as well as between the land cover characteristic driving factors. Among these factors, the interaction between FVC and ISR is the strongest. The interactions between river distance and other driving factors, as well as the interactions between slope and other driving factors, are strong in Zhuhai. Among them, the interaction between river distance and slope is the strongest. The interactions between topographic factors and water network factors are strong in Hong Kong, and the interactions between FVC and the factors of these two classes are also strong. Among them, the interaction between river density and slope is the strongest.

4.3. Spatial Regression Analysis of the Driving Mechanisms of Driving Factors on Rainstorm Waterlogging

4.3.1. The Spatial Non-Stationarity Between Rainstorm Waterlogging and Driving Factors

GWR was used to reflect and visualize the spatial differences of the relationships between different driving factors and the rainstorm waterlogging. The results are shown in Figure 7. For the land cover characteristic factors, extremely high regression coefficient values were determined for ISR, mainly distributed in the north of Dongguan, northeast of Foshan, southwest of Guangzhou, northeast of Jiangmen, middle of Shenzhen, and south of Zhongshan, occupying 2.11% of the GBA. The extremely low regression coefficient values for ISR are always distributed or staggered around the areas with extremely high regression coefficient values; the low values occupy 3.69% of the GBA. Most of the remaining areas consist of mainly negative correlations and insignificant correlations. The areas with a significant relationship between Shape_MN and rainstorm waterlogging are scattered in various cities in the GBA. The extremely high regression coefficient values for Shape_MN only occupy 0.25% of the GBA, and are mainly located in downtown Guangzhou and north Shenzhen. The extremely low regression coefficient values for Shape_MN are either staggered around the areas with extremely high regression coefficient values or scattered in various cities in the GBA, occupying in total 4.80% of the whole area. In addition, similarly to ISR, most of the areas have an insignificant negative correlation between Shape_MN and rainstorm waterlogging. In contrast from the correlation coefficients of ISR and Shape_MN, the regression coefficient for AI and rainstorm waterlogging density shows an obvious spatial pattern with a multi-core and gradually decreasing trend. The areas with extremely high regression coefficient values for AI are mainly distributed in the north of Dongguan, northeast of Foshan, southwest of Guangzhou, northeast of Jiangmen, and middle of Shenzhen, occupying 3.59% of the whole area. The regression coefficient values decrease level-by-level around the core. Furthermore, the areas with extremely low regression coefficient values are very large, occupying 8.12% of the whole area. For FVC, the areas with extremely high and low regression coefficient values are distributed in a staggered pattern and become clusters that are surrounded by large areas with no significant correlations. These clusters are mainly located in the north of Dongguan, the junction area of Foshan and Guangzhou, south Huizhou, the junction area of Jiangmen and Zhongshan, the middle of Shenzhen, the junction area of Huizhou and Shenzhen, and the junction area of Zhongshan and Zhuhai. The areas with extremely high and low regression coefficient values occupy 2.14% and 2.23% of the whole area, respectively.

For topographic factors, the extremely high regression coefficient values of elevation and slope are both mainly distributed in the north of Dongguan, northeast of Foshan, southwest of Guangzhou, and northeast of Shenzhen, occupying 0.80% and 0.92% of the whole area, respectively. Moreover, areas with positive high regression coefficient values are often surrounded by areas with negative low regression coefficient values. Topographic factors play a dominant role in the areas with negative low regression coefficient values. The extremely low regression coefficient values of elevation and slope occupy 2.77% and 4.31% of the entire area, respectively.

In most areas of the GBA, there is no significant correlation between the water network factors and rainstorm waterlogging disasters. The areas with extremely high and low regression coefficient values for the impacts of river density and river distance on rainstorm waterlogging are mainly located in the central areas of the various cities in the GBA. For river density, the extremely high regression coefficient values are mainly distributed in the south of Dongguan, northeast of Foshan, east of Guangzhou, northeast of Jiangmen, and in the middle of Shenzhen, occupying 2.81% of the whole area. The areas with extremely low regression coefficient values are either surrounded by areas with extremely high regression coefficients or distributed around the areas with extremely high regression coefficients, occupying 2.76% of the whole area. The spatial distributions of the extremely high and low regression coefficient values of the influence of the river distance on rainstorm waterlogging are basically opposite of those of the river density, accounting for 1.28% and 2.76% of the entire GBA, respectively.

4.3.2. Risk Assessment of Rainstorm Waterlogging in the GBA

Based on the analysis of the driving mechanism on rainstorm waterlogging, we found that there is collinearity between the various driving factors, which may affect the prediction accuracy of GWR model. In addition, a large number of driving factors will seriously affect the calculation efficiency of the GWR model. To overcome these limitations, before running the GWR model, we did principal component analysis on eight driving factors, and extracted three principal components (i.e., Principal Component 1, Principal Component 2, and Principal Component 3), that are unrelated with each other. The total variance explained is shown in

Table 4. The result of PCA.

Index	Principal Component 1	Principal Component 2	Principal Component 3
Total	3.810	1.185	1.003
Variance Percent	47.626	14.817	12.538
Cumulative Percent	47.626	62.443	74.981

. The cumulative percent of these three principal components is 74.981%. Moreover, for Principal Component 1, ISR, Shape_MN, and AI have high share explaining the rainstorm waterlogging degree, which reveals that this principal component mainly represents land cover characteristics. For Principal Component 2, elevation and slope have high share explaining the rainstorm waterlogging degree, which reveals that this principal component mainly represents topographic characteristics. For Principal Component 3, elevation, slope, and river distance have a high share explaining the rainstorm waterlogging degree, which reveals that this principal component mainly represents the characteristics of topography and water network.

These principal components were used to run the GWR model to evaluate the risk of rainstorm waterlogging in the GBA. Furthermore, because forestland and water areas are permeable and will not be affected by rainstorm waterlogging, these areas were masked when running the GWR model. The equal interval method was used to classify the evaluation results of rainstorm waterlogging into five categories: very low risk, low risk, medium risk, high risk, and very high risk. The estimated standard deviation of the residuals of the GWR model is 0.006. The R² and the adjusted R² are 0.9869 and 0.9849, respectively, and the difference between them is 0.002, indicating that the imitative effect of the model is credible. The range of regression coefficient of Principal component1, Principal component 2, and Principal component 3 are −0.095–0.179, −0.041–0.000, and −0.045–0.084, respectively. Therefore, Principal component 1, which represents the land cover characteristics, had the greatest contribution to the risk assessment.

Figure 8 shows the results of the risk assessment of rainstorm waterlogging in the GBA. It can be seen from the figure that the risk of rainstorms and waterlogging in the GBA presents a multi-core trend and a pattern of surrounding distribution level by level. The very high rainstorm waterlogging disaster risk areas are mainly distributed in the southwest of Guangzhou and in the central area of Shenzhen, forming two major core areas for rainstorm waterlogging disasters, which together occupy 1.67% of the entire area. High rainstorm waterlogging disaster risk areas are distributed around the two main core areas and extend to the northeast of Foshan and the northern part of Hong Kong. Three sub-core areas with high rainstorm waterlogging risk formed in the northeast of Jiangmen, east of Zhuhai, and northern Dongguan. The area with a high rainstorm waterlogging disaster risk accounts for 4.18% of the entire area. In addition to the medium-risk areas surrounding the extremely-high-risk and the high-risk areas, there are also medium-risk areas in the southeast of Zhaoqing, the north of Foshan, the middle of Huizhou, and the middle of Zhuhai. The medium-risk area occupies 8.24% of the entire area. The risk of rainstorm waterlogging in other regions is relatively low. The risk assessment results are quite similar with the actual spatial distribution of rainstorm waterlogging, which is shown in Figure 2. Therefore, the PCA-GWR can be applied to other areas in the future to predict rainstorm waterlogging risks on the basis of driving factors.

5. Discussion

5.1. Driving Mechanisms of Rainstorm Waterlogging in the GBA

The driving mechanisms of rainstorm waterlogging are different between the whole GBA and the cities in the GBA. For the whole GBA and the five cities of Foshan, Guangzhou, Huizhou, Jiangmen, and Zhaoqing located in the northern hilly area, land cover characteristic factors are the factors that are most relevant to the distribution of rainstorm waterlogging spots. Of the land cover characteristic factors, ISR, Shape_MN, and AI are positively correlated with the density of rainstorm waterlogging spots, while FVC is negatively correlated with the density. This is consistent with the research of Zhang et al. [83] and Guo et al. [84], who believe that changes in the area, shape, and fragmentation of impervious surfaces will affect the water permeability of the surface, thereby affecting the occurrence of rainstorm waterlogging. The reasonable distribution of green space in these areas also has a positive effect on reducing surface rainfall runoff. In Dongguan and Shenzhen, the developed coastal cities on the east bank of the Pearl River, the topographic factors are the strongest factors affecting rainstorm waterlogging, and both elevation and slope have negative correlations with rainstorm waterlogging density. This result is consistent with the study of Tehrany et al. [31], who concluded that high-altitude areas are less likely to experience rainstorm waterlogging events than relatively flat areas. However, for Zhongshan and Zhuhai on the west bank of the Pearl River, the water network factors have high degrees of impact on rainstorm waterlogging. The correlation between river density and rainstorm waterlogging is not significant, while river distance and rainstorm waterlogging have a relatively strong negative correlation. Wang [85] found that water network factors are negatively correlated with rainstorm waterlogging in Shanghai, but, unlike our results for Zhongshan and Zhuhai, her study revealed that the influence of water network factors on rainstorm waterlogging is lower than the influence of elevation. This difference may be due to the different levels of urban development. Shenzhen, Dongguan, and Shanghai are relatively developed, with impervious surfaces evenly distributed in flat areas. In Zhongshan and Zhuhai, construction land has mainly increased, there are dense water networks, and the urban drainage capacities have not been correspondingly improved. This combination of factors causes the pattern that the closer an area is to the water network, the denser the distribution of construction land is and the higher the risk of rainstorm waterlogging is. There is no obvious rule in the impact of driving factors on rainstorm waterlogging in Hong Kong. Hong Kong is one of the cities with the highest rainfall in the surrounding Pacific Ocean region. In addition, Hong Kong is mountainous, and its downtown area is located in a low-lying area, resulting in an extremely high risk of rainstorm waterlogging. However, rainstorm waterlogging events rarely occur in Hong Kong due to a series of flood prevention measures taken by the government [86]. Although Macao has a high urbanization rate, low terrain undulations, and abundant rainfall, it is mainly composed of islands with relatively small areas, resulting in relatively strong drainage capacity. Therefore, there are no rainstorm waterlogging events in Macao.

Simultaneously, all driving factors have enhanced the degree of the explanatory power of the spatial distribution of rainstorm waterlogging through their interactions, indicating that the frequent occurrence of rainstorm waterlogging in the region is caused by the interaction of multiple driving factors. Binary enhancement mainly occurs within the classes of driving factors, while nonlinear enhancement occurs between different classes of driving factors because the driving factors within the classes are often correlated. The interactive pairs with a high degree of explanatory power mainly include the driving factors within the land cover characteristic category and the land cover characteristic factors paired with the topographic factors. ISR and other landscape patterns of impervious surfaces have interactive influences on the occurrence of rainstorm waterlogging. The interactions between ISR and topographic factors have also been shown to enhance the power of determinants in influencing the occurrence of rainstorm waterlogging. The above interaction relationship can be described as a large impervious surface area with an irregular distribution shape, a large impervious surface area with a concentrated distribution, and a large impervious surface area with a gentle terrain slope. This shows that if an impervious surface distribution is too concentrated, it will result in the reduction of the rainwater seepage space. Moreover, the areas where rainstorm waterlogging is frequent include the older urban areas that experienced earlier development, in which the shape of the impervious surface distribution is more complicated and the terrain is relatively flat, causing the rainwater overflow to be more complicated. Thus, the rainwater fluidity is poor, so rainwater easily gathers in a short time to cause waterlogging.

5.2. Suggestions to Mitigate the Risk of Rainstorm Waterlogging in the GBA

Most of the GBA region has a subtropical monsoon climate, with a mild climate and abundant rainfall. The GBA often experiences typhoons from June to October each year; typhoons are associated with heavy rainfall, resulting in the occurrence of rainstorm waterlogging. Previous studies have confirmed that optimizing the spatial patterns of impervious surfaces can effectively reduce the surface runoff coefficient and decrease the risk of rainstorm waterlogging [17]. Increasing urban green space can also reduce surface runoff and slow rainstorm waterlogging [87]. In addition, reasonable increases in water storage infrastructure can effectively reduce the risk of rainstorm waterlogging. Based on the GD and GWR models, this paper explains the causes of rainstorm waterlogging and their spatial heterogeneity, evaluates the risk of rainstorm waterlogging, and finally gives suggestions for controlling rainstorm waterlogging risk in the study area.

This study reveals that downtown Guangzhou and the central area of Shenzhen are areas with extremely high risks of rainstorm waterlogging, although their major factors that affected the rainstorm waterlogging spot distribution are different. The contributions of land cover characteristics are more than other factors in Guangzhou, while in Shenzhen, the topographic factors contributed more than land cover characteristic factors, except FVC. The urbanization process of these areas was carried out earlier than other cities; hence, the spatial patterns of impervious surfaces are complicated, and the drainage infrastructures are outdated, which also increases the risk of rainstorm waterlogging. Quan et al. [88] and Gaitan and Veldhuis [14] believed that areas with high levels of urbanization and outdated drainage infrastructure have a higher risk of rainstorm waterlogging, and the impacts of rainstorm waterlogging disasters in these areas are also more severe. As the spatial patterns of these regions are relatively complex, with relatively large impervious surface areas, it is possible to increase the green areas and lay permeable materials to improve the rainwater absorption capacities in these areas. In addition, drainage pipe networks can be upgraded to improve rainwater drainage efficiency in areas with low terrain, where it is easy to collect rainwater.

Areas with a high risk of rainstorm waterlogging are located in northern Dongguan, northeastern Foshan, the periphery of downtown Guangzhou, northern Hong Kong, northeastern Jiangmen, most of Shenzhen, and eastern Zhuhai, with the common factors being that the surfaces have poor water permeability and the impervious surfaces are tightly distributed. Increasing green infrastructure and arranging spatial patterns of impervious surfaces in a rational manner would provide useful ways to reduce surface runoff and mitigate waterlogging in these areas. In addition, some areas are affected by relatively flat terrain and poor drainage capacity; these areas should be considered in the construction of drainage infrastructure.

For the remaining areas with low and medium risks of rainstorm waterlogging, the driving factors show weak or even no correlations with rainstorm waterlogging, except the AI of impervious surfaces, which shows a very strong negative correlation with rainstorm waterlogging. This is because, under the rapid urbanization process, the impervious surface expanded from the center to the periphery in these areas, and most of these areas are hilly landforms. Therefore, the impervious surface distribution is relatively regular, which will be helpful for water draining. Even if the impervious surface is densely distributed in some areas, the rainwater drainage efficiency is improved due to the steep terrain. As the main areas of the future urbanization process, these areas should further rationally plan the spatial structure of impervious surfaces and optimize the underground drainage pipe network system to alleviate or even effectively eliminate the regional drainage pressure and rainstorm waterlogging disasters caused by rapid urban expansion.

6. Conclusions

In this paper, we analyzed the spatial distribution characteristics of rainstorm waterlogging spots in the GBA. Then, this paper took into account the complex geographic spatial distribution characteristics of large-scale regions and analyzed the impact mechanisms of rainstorm waterlogging in GBA comprehensively. Finally, the risk of rainstorm waterlogging disasters in the GBA was assessed.

It was found that the rainstorm waterlogging spots in the GBA have strong spatial agglomeration, and are mainly concentrated in the core development areas of the GBA. The occurrence of rainstorm waterlogging in the GBA is determined by multiple driving factors. Among these factors, the land cover characteristic factors that represent the permeability of the surface are the most important factors. The rainstorm waterlogging is positively related to ISR, Shape_MN, and AI, while negatively related to FVC, elevation, slope, and river distance. The correlation between river density and rainstorm waterlogging is not strong. Moreover, all of the driving factors enhance the degree of rainstorm waterlogging impacts through their interactions, and most of the interactions include the land cover characteristic factors. In addition, the impacts of driving factors on rainstorm waterlogging are different in different spatial positions, which reflects the spatial non-stationarity of the relationship between rainstorm waterlogging and driving factors. Based on the risk assessment of rainstorm waterlogging by integrating multiple driving factors, we find that the distribution of the risk of rainstorm waterlogging in the GBA presents multi-core and layer-by-layer trends. Furthermore, the southwest of Guangzhou and the central area of Shenzhen have the highest risks of rainstorm waterlogging.

These findings provide a reference for the prevention and control of rainstorm waterlogging in GBA. Meanwhile, it provides a foundation for the further research of rainstorm waterlogging in urban agglomeration. However, there are still the following problems. In this paper, we mainly consider the impact of land cover, topography, and water system distribution on rainstorm waterlogging. The inclusion of more data, such as rainfall and underground drainage data, would help to comprehensively assess the risk of rainstorm waterlogging in the GBA.