Identification of Key Determinants Influencing the Spatiotemporal Heterogeneity of Urban Resilience

1. Introduction

Urban resilience (UR) refers to the capacity of urban systems, including their social technology and social ecological networks, to rapidly recover from disruptions that occur on a specific spatiotemporal scale (Hamilton, 2009; Jabareen, 2013; Burton, 2015). It encompasses the ability to recover from disasters and other crises, while continuing to fulfill industrial, commercial, residential, governmental, and social aggregation functions (Kevin, & Trevor, 2013; Meerow, Newell, & Stults, 2016). The heterogeneity of resilience across different cities arises due to variations in natural, economic, social, and ecological factors (Adger, 2000). Investigating this heterogeneity, as well as its key determinants and the interactions among factors within and between subregions, could determine the mechanisms underlying UR. The results could then be used to inform urban planning and safety governance.

The concept of resilient cities was initially proposed by the International Council for Local Environmental Initiatives in 2002 (Cutter, Barnes, Berry, Burton, Evans, et al. 2008). The Rockefeller Foundation’s project, "100 Resilient Cities", launched in 2013, aimed to encourage governments worldwide to develop resilient urban areas (Marjolein, & Bas, 2017; Daniela, Mauricio, Roberto, & Cristóbal, 2018). In accordance with the 2030 Agenda for Sustainable Development, the construction of inclusive, safe, resilient, and sustainable cities is imperative (Hamilton, 2009). Highly resilient cities are characterized by diversity, modularity, innovation, adaptability to change, rapid feedback ability, social capital reserves, and ecosystem service capabilities (Allan, 2011). Urban elements with high resilience should possess multifunctionality and redundancy, while incorporating ecological and social diversity (Adger, Hughes, Folke, Carpenter, & Rockstrom, 2005). Additional important factors include multi-scale network connectivity, as well as adaptive planning and construction strategies (Ahern, 2011).

The research methods used in studies of UR primarily include function modeling, threshold analysis, social network analysis, scenario analysis, layer overlay, and the geographical detector (geodetector) (Bruneau, Chang, Eguchi, Lee, O’Rourke, et al. 2012; Wang, & Xu, 2017). Function modeling is a straightforward and user-friendly approach (Du, Zhang, Wang, Tao, Li., & Li, 2020); however, it relies on limited datasets and few parameters, leading to significant evaluation errors (Raquel, Leire, & Josune, 2017). The threshold method enables the analysis of complex system evolution and coordination, but is sensitive to threshold selection and noise, resulting in coarse segmentation outcomes (Kamila, Peter, & Porfirio, 2018). A social network analysis can examine the connectivity of social networks and the structural resilience between nodes (Christmann, & Ibert, 2012); however, it overlooks isolated nodes as well as subjective motivations, functions, and the nature of social networks. A scenario analysis enables the impact of multiple factors on UR to be simulated under different scenarios (Josune, Patricia, Raquel, Jose, & Leire, 2019); however, scenario assumptions often rely on personal subjective judgments and biases that reduce prediction accuracy. The layer overlay method intuitively expresses spatial differences in UR by overlaying multiple element layers within a geographic information system (GIS) environment (Meerow, & Newell, 2016); however, identifying interaction effects among multiple factors has proven challenging. The geodetector statistical tool quantitatively analyzes the mechanisms influencing each driving factor, as well as the interactions between factors (Wang, Li, Christakos, Liao, Zhang, Gu, et al., 2010), although its spatial expression capability remains weak (Zhang, Dou, Liu, & Gong, 2023).

The study of UR has evolved from a conceptual framework to a measurable construct, encompassing multiple disturbances rather than singular events. This has shifted the focus of study from infrastructure-centric construction to the humanistic aspects of resilience, and has involved a transitioning from static notions of resilience to dynamic processes (Walker, Holling, Carpener, & Kinzig, 2004; Carvalho, Martins, Marta-Almeida, Rocha, & Borrego, 2017). However, investigations into the spatiotemporal heterogeneity of UR often concentrate solely on spatial variations across the entire study area (Li, Kappas, & Li, 2018). They have rarely considered the interaction of multiple factors or the identification of key determinants. Additionally, the interaction among multiple factors; identification of key determinants; and the contributions made by intra- and inter- subregions toward the spatiotemporal heterogeneity of UR have been neglected. Therefore, it is imperative to fully determine the mechanisms underlying the spatiotemporal heterogeneity of UR.

A UR index system and evaluation model were constructed in this study. A GIS was used to visually analyze the spatiotemporal evolution and driving factors of UR. The analysis of cold and hot spots revealed the agglomeration characteristics of UR. The spatial and dynamic evolution of UR was examined by analyzing its spatial directionality using a standard deviation ellipse. Furthermore, the contributions of spatiotemporal heterogeneity within and between subregions were analyzed using the Theil index. The gedetector tool was used to identify how the key determinants of UR and their interactions affected its spatiotemporal heterogeneity. This study provides a theoretical basis and practical support for evaluating, promoting, and planning UR. Noteworthy innovations of the study include: (1) proposing a factor interaction increment; (2) conducting an in-depth analysis on factors contributing to the spatiotemporal heterogeneity of UR; (3) systematically expressing various aspects of UR, such as its dynamic spatiotemporal evolution, key determinants, agglomeration characteristics, contributions from intra- and inter-subregions, and single- and two-factor interactions.

2. Materials and Methods

2.1. Study Area

Hunan Province, situated in the middle reaches of the Yangtze River in China, covers a total land area of 211,829 km² between 108°47′ E and 114°15′ E longitude and 24°38′ N and 30°08′ N latitude. The province has a humid monsoon climate within the middle subtropical zone, characterized by an average annual temperature of 17.7°C and an average annual rainfall of 1200−1700 mm. Its altitude above sea level ranges from -123 to 2093 m. Most of the mountains with altitudes exceeding 1000 m are located in eastern, southern, and western Hunan Province. The central region is marked by undulating hills interspersed with valley basins at altitudes below 500 m. The Dongting Lake Plain is located in northern Hunan, where altitudes are mostly below 50 m (Figure 1). Numerous rivers flow through Hunan Province, including Xiangjiang River, Zijiang River, Yuanjiang River, and Lishui River, which drain into Dongting Lake and Yangtze River from south to north.

Hunan Province is traditionally divided into four subregions based on the spatial distribution of its cities and the layout of economic and social development. The northern region encompasses Yueyang, Changde, and Yiyang. The central region includes Changsha, Zhuzhou, Xiangtan, Shaoyang, and Loudi. The southern region comprises Hengyang, Chenzhou, and Yongzhou. Finally, the western region consists of Huaihua, Zhangjiajie, and Xiangxi.

The cities in the four traditional subregions each possess unique advantages in terms of their location, resources, and industries. However, their overall resilience is hindered by environmental degradation, resource depletion, industrial homogenization, uneven regional development, and spatial conflicts between production, livelihoods, and ecology. Identifying the key determinants of the spatiotemporal heterogeneity of UR would provide theoretical support for urban planning and resilience enhancement within the Hunan and Yangtze River economic belt.

2.2. Construction of an UR Index System

Nine indexes were established to assess the economic, social, and ecological resilience, considering the relative and absolute differences in the urban development of Hunan Province. These indexes provided a comprehensive framework to gauge UR across the province (Table 1).

The selection basis of each index was as follows. Per capita gross domestic product (GDP) was used to reflect the overall economic strength and pressure resistance of a city. The proportion of tertiary industry was used to assess the influence of the rationality and diversity of the urban economic structure on various factors. Per capita social consumption was adopted to reflect the market scale and economic development of a city. The average wage of on-the-job employees was used to measure the overall living standard of urban residents. The number of hospital beds per 10,000 people and doctor density were used to reflect the adaptability and recovery capacity of the medical response to disasters in urban areas. Per capita green land in park areas and the greening coverage rate of the built-up area were adopted to reflect the ability of a city to restore and reconstruct the urban ecosystem in response to ecological challenges and disturbances. The sewage treatment rate was used to measure the capacity of urban areas to handle pollutants and protect the environment.

The values of the UR evaluation indexes were mainly obtained from the China Statistical Yearbook (2015−2023), the Hunan Statistical Yearbook (2015−2023), the statistical yearbooks of cities, and the social and economic development bulletins of cities in Hunan Province. Panel data for 2014, 2018, and 2022 were used to analyze the spatiotemporal heterogeneity of UR.

2.3. Construction of the UR Evaluation Model

The entropy weight method was adopted in this study to determine the weight of each index. It is an objective weighting approach with a simple operation and high degree of reliability (Sherrieb, Norris, & Galea, 2010). An index has a large degree of variation, contains a large amount of information, and has a large weight if the information entropy is small.

(1) Standardization of evaluation indexes

Suppose there are n cities and m evaluation indexes for each city, then matrix X is constructed as follows:

$$X=\left[\begin{array}{cccc}{x}_{11}& x_{12}& \cdots & x_{1m}\\ x_{21}& x_{22}& \cdots & x_{2m}\\ ⋮& ⋮& \ddots & ⋮\\ x_{n1}& x_{n2}& \cdots & x_{nm}\end{array}\right],$$

(1)

y_ij (i = 1, 2, …, n; j = 1, 2, …, m) is obtained after the standardization of the x_ij because the indexes used in this study were all positive.

$$y_{ij}={\displaystyle \frac{x_{ij}-\min \left(x_{1j},x_{2j},\cdots x_{nj}\right)}{\max \left(x_{1j},x_{2j},\cdots x_{nj}\right)-\min \left(x_{1j},x_{2j},\cdots x_{nj}\right)}},$$

(2)

(2) Calculation of index proportion

The proportion p_ij is calculated according to the ratio of y_ij to the sum of y_ij (i = 1, 2, …, n; j = 1, 2,…, m).

$$p_{ij}={\displaystyle \frac{y_{ij}}{{\displaystyle {\sum }_{i=1}^n{y}_{ij}}}},$$

(3)

(3) Calculation of the information entropy e_j of the jth (j = 1, 2,…, m) index.

$$e_j=-{\displaystyle \frac{1}{\ln n}}{\displaystyle \sum _{i=1}^n{p}_{ij}}\ln p_{ij},$$

(4)

(4) Calculation of the weight w_j of the jth (j = 1, 2,…, m) index.

$$w_j={\displaystyle \frac{1-e_j}{{\displaystyle {\sum }_{j=1}^m\left(1-e_j\right)}}},$$

(5)

(5) The UR evaluation model was constructed based on a comprehensive index method to calculate the resilience UR_i of the ith city.

$$U{R}_i={\displaystyle \sum _{j=1}^m{w}_j{y}_{ij}},$$

(6)

(6) Hierarchical classification of UR

The UR was classified into five grades in accordance with the distribution characteristics of UR in Hunan Province, namely, the lowest resilience if UR ∈ [0, 0.15), low resilience if UR ∈ [0.15, 0.2), moderate resilience if UR ∈ [0.2, 0.3), high resilience if UR ∈ [0.3, 0.4), and the highest resilience if UR ∈ [0.4, 1.0].

2.4. Cold and Hot Spots of UR

Hot spot analysis can identify the spatial location of the clustering of high or low values of UR based on z scores and p values. The p value represents the probability that the observed spatial pattern is constructed by a random process. The z-score is a multiple of a standard deviation. For instance, the z-score is 2.7 times the standard deviation if the z-score = +2.7. The z-score is calculated by G_i* as follows:

$$G_i^{*}={\displaystyle \frac{{\displaystyle \sum _{j=1}^n{w}_{ij}{x}_j}-\overline{X}{\displaystyle \sum _{j=1}^n{w}_{ij}}}{S\sqrt{\left[n{\displaystyle \sum _{j=1}^n{w}_{ij}^2}-{\left({\displaystyle \sum _{j=1}^n{w}_{ij}}\right)}^2\right]/\left(n-1\right)}}},$$

(7)

$$\overline{X}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{j=1}^n{x}_j},$$

(8)

$$S=\sqrt{{\displaystyle \frac{1}{n}}{\displaystyle \sum _{j=1}^n{x}_j^2}-{\left(\overline{x}\right)}^2},$$

(9)

where i is the ith city, j is the jth city in the neighborhood of the ith city, x_j is the resilience of the jth city in the neighborhood, w_ij is the spatial distance between the ith and jth cities, n is the number of cities in the neighborhood, $\overline{x}$is the average value of x_j (j = 1, 2,…, n), and S is the standard deviation.

The higher the positive z score with statistical significance, the closer the clustering of high values (hotspots). The lower the negative z score with statistical significance, the closer the clustering of low values (cold spots).

2.5. Theil Index

The Theil index was used to measure the differences in UR in the entire study area, and between and within the subregions of the study area. The greater the value of Theil index, the larger the spatial difference of UR.

$$T={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n\frac{x_i}{\overline{x}}}\ln \left({\displaystyle \frac{x_i}{\overline{x}}}\right),$$

(10)

where T ∈ [0, 1] is the Theil index values used to describe the total difference of UR, n is the number of cities, x_i is the resilience of the ith city, and $\overline{x}$is the average value of the resilience of all cities.

The T value was further divided into T_w to express the intraregional difference of UR and T_b to indicate the interregional difference of UR.

$$T=T_w+T_b,$$

(11)

$$T_w={\displaystyle \sum _{p=1}^m\left(\frac{n_p}{n}\frac{{\overline{x}}_p}{\overline{x}}\right)}{T}_p,$$

(12)

$$T_b={\displaystyle \sum _{p=1}^m\left(\frac{n_p}{n}\frac{{\overline{x}}_p}{\overline{x}}\right)}\ln \left({\displaystyle \frac{{\overline{x}}_p}{\overline{x}}}\right),$$

(13)

where T_w is the spatial heterogeneity of UR within a subregion of the study area, T_b is the spatial heterogeneity of UR between subregions, m is the number of subregions of the study area, n_p is the number of cities in the pth (p = 1, 2, …, m) subregion, ${\overline{x}}_p/\overline{x}$ is the ratio of the average of UR in the pth subregion to the average of UR of all cities in Hunan Province, and T_p is the Theil index of the pth subregion.

2.6. Factor Detection of UR

Geodetector was used in this study to analyze the factors influencing the spatiotemporal heterogeneity of UR, due to its ability to detect the spatial differentiation, driving forces, and interactions of geographical phenomena (Wang, & Hu, 2012). The k-means clustering method was initially applied to convert numerical variables into discrete ones, i.e., the domain of variable x was partitioned into h layers (Zhu, Yu, Chen, Tan, & Yuan, 2024). Subsequently, the q value was adopted to detect the spatial heterogeneity of y, i.e., the extent to which the detection factor x affects the spatial differentiation of y. The value of q was calculated as follows:

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

(14)

where q∈ [0, 1] represents the degree of influence of the detection factor x on UR at the 95% confidence level, h denotes the hth layer (h = 1, 2, …, L) where the value of factor x is equally divided into L intervals, N and N_h signify the number of samples in the study area and the detection area respectively, σ² and σ_h² indicate the variances of y in the study area and the detection area respectively, and ${\displaystyle \sum _{h=1}^L{N}_h{\sigma }_h^2}$is the sum of variances within the hth layer.

A straightforward transformation of the q value conformed to the noncentral F distribution as follows:

$$F={\displaystyle \frac{N-L}{L-1}}{\displaystyle \frac{q}{1-q}}~F\left(L-1,N-L;\lambda \right),$$

(15)

$$\lambda ={\displaystyle \frac{1}{{\sigma }^2}}\left[{\displaystyle \sum _{h=1}^L{\overline{Y}}_h^2}-{\displaystyle \frac{1}{N}}{\left({\displaystyle \sum _{h=1}^L\sqrt{N_h}{\overline{Y}}_h}\right)}^2\right],$$

(16)

where λ is a noncentral parameter, and ${\overline{Y}}_h$is the average value of UR in the hth layer.

The larger the q, the stronger the degree of influence of factor x on y; conversely, the smaller the q, the weaker the influence. The factor x fully governs the spatial distribution of y when q = 1. In contrast, factor x has no association with y when q = 0.

2.7. Detection of Interactions among UR Factors

The detection of interactions among UR factors can identify whether the influence of an interaction between the factors x₁ and x₂ will enhance or diminish the degree of influence on the dependent variable y, or whether the influence of these factors on y is independent. The q values of the factors x₁ and x₂ for y, namely q(x₁) and q(x₂), respectively, were calculated. The q value for the interaction between x₁ and x₂ was then calculated, i.e., q(x₁∩x₂) was computed, which is a new polygon intersected by two vector layers of x₁ and x₂. Finally, q(x₁), q(x₂), and q(x₁∩x₂) were compared based on formula (17) to detect the type of interaction between the two factors (T_id(q)).

$$T_{id}\left(q\right)=\left\{\begin{array}{ll}\mathrm{Nonlinear}\mathrm{weakening}\hfill & \mathrm{if}q\left(x1\cap x2\right)<\mathrm{Min}\left(q\left(x1\right),q\left(x2\right)\right)\hfill \\ \mathrm{Sin}\mathrm{gle}-\mathrm{factor}\mathrm{nonlinear}\mathrm{weakening}\hfill & \mathrm{if}\mathrm{Min}\left(q\left(x1\right),q\left(x12\right)\right)<q\left(x1\cap x2\right)<\mathrm{Max}\left(q\left(x1\right),q\left(x2\right)\right)\hfill \\ \mathrm{Two}-\mathrm{factor}\mathrm{enhancement}\hfill & \mathrm{if}q\left(x1\cap x2\right)>\mathrm{Max}\left(q\left(x1\right),q\left(x2\right)\right)\hfill \\ \mathrm{Independence}\hfill & \mathrm{if}q\left(x1\cap x2\right)=q\left(x1\right)+q\left(x2\right)\hfill \\ \mathrm{Nonlinear}\mathrm{enhancement}\hfill & \mathrm{if}q\left(x1\cap x2\right)>q\left(x1\right)+q\left(x2\right)\hfill \end{array}\right.,$$

(17)

Nonlinear weakening in formula (17) implies that the interaction of two factors nonlinearly reduces the influence of variables. Single-factor nonlinear weakening indicates that an interaction between two factors weakens the influence of variables. Nonlinear enhancement implies that the degree of influence of two interacted factors is greater than the sum of the degrees of influence of the two factors acting independently. Two-factor enhancement means that the degree of influence of two interacted factors is greater than that of each factor acting independently, and its influence is less than that of nonlinear enhancement.

Factor interaction increment (I_i(q)) is proposed in this study, which is defined as the difference in the degree of influence produced by the interaction between two factors minus that produced by one factor. Factor interaction increment can effectively analyze the contribution of factor interaction on the improvement of UR, and it is calculated as follows:

$$I_i\left(q\right)=\left\{\begin{array}{ll}-\left(\mathrm{Min}\left(q\left(x1\right),q\left(x2\right)\right)-q\left(x1\cap x2\right)\right)\hfill & \mathrm{if}\mathrm{nonlinear}\mathrm{weakening}\hfill \\ -\left(\mathrm{Max}\left(q\left(x1\right),q\left(x2\right)\right)-q\left(x1\cap x2\right)\right)\hfill & \mathrm{if}\sin \mathrm{gle}\mathrm{factor}\mathrm{nonlinear}\mathrm{weakening}\hfill \\ q\left(x1\cap x2\right)-\mathrm{Max}\left(q\left(x1\right),q\left(x2\right)\right)\hfill & \mathrm{if}\mathrm{two}-\mathrm{factor}\mathrm{enhancement}\hfill \\ 0\hfill & \mathrm{if}\mathrm{independence}\hfill \\ q\left(x1\cap x2\right)-q\left(x1\right)-q\left(x2\right)\hfill & \mathrm{if}\mathrm{nonlinear}\mathrm{enhancement}\hfill \end{array}\right.$$

(18)

3. Results and Analysis

3.1. Temporal Evolution of UR

The weights of the UR evaluation indexes in Hunan Province were determined and shown in Table 1. The resilience levels of 14 cities in Hunan Province in 2014, 2018, and 2022 were calculated according to the UR evaluation model (Figure 2). The average values of UR in Hunan Province were 0.2692 in 2014, 0.3021 in 2018, and 0.3422 in 2022, which represented an increase of 27% from 2014 to 2022, demonstrating an ascending trend. This was because the economic and social development, infrastructure, public services, industrial transformation, and innovation capabilities of cities in Hunan Province were constantly improving over the study period.

The resilience of 11 cities increased from 2014 to 2022, with the cities experiencing a significant increase including Xiangtan, Hengyang, Yueyang, and Huaihua. In recent years, China has implemented various policies and national strategies, such as the strategy to enhance the central China region, which has been achieved through increasing infrastructure construction, ensuring and improving people's livelihoods, and building an ecological civilization. Additionally, investment into Hunan Province in terms of the urban economy, society, and ecology has continuously increased, and the levels of urban safety, health, and sustainable development have risen year by year. Therefore, the UR of most cities in Hunan Province has also increased year by year.

Changsha, Zhuzhou, and Xiangtan maintained a high UR from 2014 to 2022. This was mainly due to the Changsha−Zhuzhou−Xiangtan economic integration that has been promoted by the Hunan provincial government since 1997. The Changsha−Zhuzhou−Xiangtan urban agglomeration was initially established, and has since become the core growth pole of economic development in Hunan Province. Changsha is the political, economic, cultural, educational, and financial center of Hunan Province, with a high level of urbanization and the highest UR.

Changsha and Yiyang have undergone a transition from a high UR in 2014, to a low UR in 2018, and the highest UR in 2022. This was because the profit growth of industrial enterprises above a designated size in Changsha declined, and the output of various heavy industry sectors decreased sharply after 2015, largely due to the relocation of Sany Heavy Industry Co., Ltd. from Changsha to Beijing. After 2018, industrial clusters were established in Changsha, including engineering machinery, food processing, automobiles, new materials, electronic information, biomedicine, and culture and tourism. As a result, Changsha's economy has experienced rapid development. Catastrophic flood disasters in Yiyang City affected 1.47 million people in 2016, with a direct economic loss of 2.6 billion CNY, and a further 0.78 million people in 2017, with a direct economic loss of 2.15 billion CNY. Consequently, Yiyang's UR declined in 2018. After 2018, industries in Yiyang, such as high-end equipment manufacturing, food processing, electronic information, and new materials, have developed rapidly. This has resulted in a rapid increase in UR from 2018 to 2022.

The UR in Changde was low in 2014, highest in 2018, and high in 2022. This is because the development of highly-polluting, high-energy consumption industries, such as coal and chemicals, has been restricted by environmental protection policies in China after 2018. The UR in Hengyang was lowest in 2014, with a significant increase in 2018, and was highest in 2022. This is because per capita GDP, per capita social consumption, per capita green area, and the sewage treatment rate have increased rapidly in recent years.

The UR in Zhangjiajie declined significantly from 2018 to 2022. This was mainly because specific measures, such as closing public places and prohibiting the movement of people due to the corona virus disease 2019 (COVID-19) from 2019 to 2022, restricted the development of the tourism industry, which is the pillar industry of Zhangjiajie.

3.2. Spatial Evolution of UR

The resilience of 14 cities in Hunan Province was spatially visualized using GIS to determine the spatial heterogeneity of UR (Figure 3).

There were obvious spatial differences in UR in Hunan Province, namely, a high resilience in the northeastern and low resilience in the southwest Hunan. In 2014, Changsha has the highest resilience, Xiangtan and Zhuzhou were high-resilience cities, while Shaoyang, Hengyang, and Loudi were the lowest-resilience cities.

In 2018, Changsha and Zhuzhou were the highest-resilience cities. The high-resilience cities included Zhangjiajie, Changde, and Xiangtan. Among them, Zhangjiajie and Changde transitioned from moderate-resilience cities in 2014 to high-resilience cities in 2018. Hengyang upgraded from a lowest-resilience city in 2014 to a moderate-resilience city in 2018 due to its rapid economic, social, and ecological development. In contrast, Yiyang was downgraded from a moderate-resilience city in 2014 to a low-resilience city in 2018. Shaoyang was always the lowest resilience city from 2014 to 2022.

In 2022, the highest-resilience cities were Changsha, Xiangtan, and Zhuzhou. The high-resilience cities included Changde, Yueyang, Yiyang, Hengyang, and Chenzhou. Huaihua and Loudi were upgraded from low-resilience cities in 2018 to moderate-resilience cities in 2022. Yiyang was upgraded from a low-resilience city in 2018 to a high-resilience city in 2022. Zhangjiajie was downgraded from a high-resilience city in 2018 to a low-resilience city in 2022, with a significant decline.

Therefore, the UR of Hunan had specific spatiotemporal heterogeneity, spatial agglomeration characteristics, and polarization phenomenon. The UR of the Changsha−Zhuzhou−Xiangtan urban agglomeration was higher than that of other cities. The highest-resilience cities had a positive spillover effect on their surrounding cities. Over time, the number of cities that upgraded from low to high resilience increased; the degree of agglomeration of high-resilience cities gradually increased; and the gap in resilience between cities narrowed.

The spatial heterogeneity of UR in Hunan Province is influenced by multiple factors. The northeastern cities, such as Changsha, Zhuzhou, and Xiangtan, had high resilience due to the unique policies that have been implemented within them, such as economic integration, and the massive investment in social public services and infrastructure. In contrast, the southwestern cities, such as Shaoyang, had low resilience because of their high level of dependence on resources, and their poor infrastructure, economic development, and social public services.

3.3. Directivity Analysis of the Spatial Distribution of UR

Ellipse parameters, such as the ellipse area, barycenter, semi-minor axis, semi-major axis, oblateness, and azimuth, were obtained according to the UR of Hunan Province in 2014, 2018, and 2022 using the standard deviation ellipse method. These parameters can be used to further explore the process and direction of UR evolution spatially and dynamically (Table 2 and Figure 4).

The areas of the ellipses, which represent the extent of the spatial distribution of UR, were 68,098.97 km² in 2014, 70,299.88 km² in 2018, and 69,548.45 km² in 2022, respectively. This indicates that the spatial impact of UR was largest in 2018, followed by 2022, and was smallest in 2014.

The standard deviation ellipse exhibited a distribution pattern from northwest to southeast in Hunan Province. The semi-major axis, which indicated the direction of the dense spatial distribution of UR, initially increased from 2014 to 2018 and subsequently decreases from 2018 to 2022. This suggests that the UR displayed a trend of initially strengthening and then weakening, as well as initially dispersing and then aggregating, in the main axis direction from the northwest to southeast of Hunan Province.

The semi-minor axis of the ellipse, which represents the direction of the sparse spatial distribution of UR, gradually increased from 2014 to 2022. This indicated that the pulling effect of UR on the standard deviation ellipse in the northeast and southwest of Hunan Province was gradually intensifying. The spatial distribution pattern of UR was relatively stable, and the spatial spillover effect of UR was not prominent.

Oblateness, which represents the directionality and centripetal force of UR, was largest (0.179) in 2018, indicating that strongest directionality of UR occurred in 2018. The oblateness of UR in Hunan Province was smallest (0.079) in 2022, suggesting the weakest directionality of UR, with the length of the semi-major axis approximating that of the semi-minor axis.

The barycenter of UR in Hunan Province was located at the junction of Changsha and Xiangtan, within 112°8'46" E to112°12'36" E and between 27°55'59" N to 27°59'25" N. The barycenter of UR migrated from the northeast in 2014 to the southwest in 2018, with a migration distance of 6.2 km. This was because in the southern Hunan region there was a vigorous promotion of industrial transformation, and upgrading, logistics, and financial expertise were obtained from the economically and commercially developed Pearl River Delta region. In the western Hunan region, an ecological and cultural tourism area was established. The barycenter of UR migrated from the northwest in 2018 to the southeast in 2022, with a migration distance of 5.50 km. This was because Hengyang and Chenzhou, situated in the southeast of Hunan Province, emerged as new growth poles for industrial transfer from the Pearl River Delta.

The azimuth of the ellipse, which represents the main trend in the direction of the spatial distribution of UR, rotated clockwise from 2014 to 2018. This indicates that the pulling effect of cities in the southwest of Hunan Province on the standard deviation ellipse gradually became stronger than that of cities in the northeast at this stage. The azimuth rotated counterclockwise from 2018 to 2022, suggesting that the ability of cities in the northeast to shape the UR spatial pattern was enhanced at this stage.

3.4. Hot Spot Analysis of UR

The distribution of hot and cold spots of UR in Hunan Province was obtained through the use of the Getis-Ord Gi* statistic in GIS (Figure 5). The UR displayed a trend of a gradual decline in hot spots and a gradual increase in cold spots from the northeast to southwest of Hunan Province.

In 2014, the hot spots of UR in Hunan Province were distributed in Changsha, Zhuzhou, and Yueyang. There were also sub-hot spots distributed in Xiangtan and Yiyang. Sub-cold spots were distributed in Xiangxi, Loudi, and Chenzhou. The cold spots were distributed in Yongzhou, Shaoyang, and Huaihua. Generally, the northeast region was dominated by hot spots, and the southwest region was dominated by cold spots.

In 2018, the hot spots of UR were distributed in Zhuzhou, located in the east of Hunan Province. The sub-hot spots were distributed in Changsha, Yiyang, Xiangtan, and Yueyang, in the northeast of Hunan Province. The sub-cold spots were located in Yongzhou and Xiangxi, and the cold spots were distributed in Shaoyang, Huaihua, and Loudi in the southwest. In general, the hot spots of UR in 2018 were mainly concentrated in the northeast of Hunan Province, and the cold spots were mainly concentrated in the southwest.

In 2022, hot spots of UR were identified in Changsha and Zhuzhou, situated in the eastern region of Hunan Province. Sub-hot spots were observed in the periphery of the hot spots, including Yueyang, Yiyang, and Xiangtan. Conversely, sub-cold spots were found in Yongzhou, Loudi, Changde, Zhangjiajie, and Xiangxi. Furthermore, cold spots were evident in Shaoyang and Huaihua in the southwest. Consequently, the distribution of cold and hot spots of UR in Hunan Province mirrored that of UR overall. This spatial pattern confirmed that northeastern Hunan Province was an agglomeration area for hot spots, while the southwest served as an agglomeration area for cold spots; thus, further validating the spatiotemporal heterogeneity of UR in Hunan Province.

3.5. The Spatiotemporal Differences in UR

The overall differences in UR in Hunan Province were measured using the Theil index (Table 3). The Theil index values of UR were 0.17026 in 2014, 0.11703 in 2018, and 0.09327 in 2022, respectively, i.e., gradual decrease. This indicates that the spatial differences in the UR in Hunan Province from 2014 to 2022 became increasingly smaller. This might be related to the strategies to enhance the central China region, the scientific approach taken toward development, and the coordination of regional development.

The spatial variances of UR between and within subregions of the traditional zones, as well as the hot spots of UR in 2014, 2018, and 2022 in Hunan Province were determined using the Theil index (Figure 6). The traditional zoning encompassed four subregions: northern, central, southern, and western Hunan. The zones based on the urban hot spots in 2014, 2018, and 2022 comprised five subregions: hot-spot, sub-hot spot, non-significant, sub-cold spot, and cold-spot areas, respectively (Table 3 and Figure 4).

The Theil index values of resilience within the subregions of Hunan (T_w), derived by the four zoning approaches, progressively diminished from 2014 to 2022. This indicated that the spatial disparity of UR within subregions was shrinking at an increasing rate, which was also the case for the entire Hunan Province. The Theil indexes of UR between subregions (T_b) obtained from the zoning approaches based on the hot spots of UR in 2014 and 2022 gradually reduced from 2014 to 2022. This indicated that the spatial variance of UR between subregions gradually narrowed. The Theil index values of UR between subregions (T_b) in 2014 were greater than those in 2018 and 2022, while the Theil index values of UR between subregions (T_b) in 2018 were nearly equivalent to those in 2022, according to the methods used for conventional zoning and the zoning of hotspots of UR in 2018. This indicated that the spatial variance of UR between subregions in 2014 was more significant than in 2018 and 2022.

There were similarities and dissimilarities in the contribution rates of the spatiotemporal variances of UR within and between subregions due to the diverse zoning methods used. Using the conventional zoning method, the Theil index values of UR within subregions were considerably larger than those between subregions. This indicated that the regional variances of UR mainly arose from the variances within the subregions. For the zoning methods based on the hotspots of UR in 2014 and 2018, the Theil indexes of UR within subregions were analogous to those between subregions. This indicated that the spatial variance of UR within subregions was approximately the same as that between subregions. Using the zoning method based on the hotspots of UR in 2022, the Theil index values of UR within subregions were significantly smaller than those between the subregions of cities. This indicated that the spatial variance in UR mainly originated from the interregional UR.

There were similarities and dissimilarities in the spatiotemporal heterogeneity of UR within the subregions obtained by the different zoning methods (Table 3 and Figure 7). Using the conventional zoning method, the variances within subregions in central Hunan in 2014, 2018, and 2022 were conspicuously higher than those in other subregions. This indicated that the UR in central Hunan exhibited the greatest spatial heterogeneity and imbalance (Figure 7a), with the UR of Changsha City in central Hunan being consistently the highest, while that of Shaoyang was always the lowest throughout the study period. Additionally, the Theil index of UR in central Hunan gradually declined from 2014 to 2022, suggesting that the spatial variance of UR in central Hunan gradually decreased. The spatial variances of UR in northern Hunan were largest in 2018 and smallest in 2022. The spatial variances of UR in western Hunan underwent negligible change from 2014 to 2022. The Theil index values of UR in southern Hunan were largest in 2014, and very small in 2018 and 2022, indicating that the spatial variance of UR in southern Hunan was substantial in 2014, and then relatively balanced in 2018 and 2022.

Using the zoning based on hotspots of UR in 2014, the internal variances in the hotspot area were significantly higher than in other subregions during the study period. This indicated that the UR of the hotspot area presented the largest spatial variance and imbalance (Figure 7b). This was because Changsha had the highest UR, while Yueyang had a low resilience, resulting in the polarization of UR during the study period. The Theil index of UR in the sub-hot spot subregion was highest in 2018, and low in 2014 and 2022. The spatial variance of UR in the non-significant subregion was largest in 2014, and small in 2018 and 2022. The spatial variance of UR in the sub-cold spot subregion was very small during the study period. The spatial variance of UR in the cold spot subregion was largest in 2022, and low in 2014 and 2018.

Using the zoning based on the hotspots of UR in 2018, the Theil index was 0 in each of 2014, 2018, and 2022 because only Zhuzhou City was located in the hotspot subregion, suggesting a lack of internal differences in UR. The Theil index values of UR in the sub-hotspot subregion were significantly higher than in other subregions during the study period. This suggests that the UR in the sub-hotspot subregion exhibited substantial spatial differences and imbalances (Figure 7c). These were attributed to the pronounced polarization between high UR areas, such as Changsha and Zhuzhou, and low UR areas, such as Yiyang during the study period. Furthermore, the spatial differences in UR were greatest within non-significant subregions in 2014 but were smallest in these areas by 2018. Similarly, the spatial difference of UR within sub-cold spot subregions peaked in 2014 but was minimized by 2022; while the differences within cold spots reached their maximum disparity by 2022, but were reduced to a minimum by 2018.

Using the zoning based on the hotspots of UR in 2022, the internal differences of UR within hotspot subregions were found to be significantly larger than those in other subregions. Moreover, minimal changes were observed in the internal differences of UR during the study period (Figure 7d). The spatiotemporal heterogeneity of UR within sub-hot spots, sub-cold spots, and cold spots exceeded those shown in Figure 6a, b, and c. The Theil index values for UR within sub-hot spots and sub-cold spots peaked in 2018 but decreased notably by 2014 and 2022. The spatial disparities of UR within cold spots reached their peak in 2022, with a minimum level in 2018. Furthermore, the spatial disparities of UR within non-significant subregions surpassed those within hotspots in 2014, but declined substantially by 2018, with the decline continuing into 2022.

In summary, among the four zoning methods, the largest Theil index was observed in the area centered on Changsha, and the spatial differences of UR decreased from 2014 to 2022. The Theil index values in other subregions were generally small.

3.6. Factor Detection Analysis

The driving factors of UR in Hunan Province were analyzed using geodetector (Figure 8). The selected factor values were categorized into five levels: highest, high, medium, low, and lowest based on the classification standard of UR. The degree of influence q and factor explanatory power p of each factor on UR were then calculated using the geodetector tool. A higher value of q indicated a larger influence on UR, while a lower value of p suggested a greater explanatory power for influencing UR.

Among the urban economic resilience indexes (x1, x2, x3), except for the proportion of tertiary industry, per capita GDP, and per capita social consumption had significantly higher values of q (> 0.9) from 2014 to 2022, making them important factors affecting the spatial heterogeneity of UR. The p value obtained for the proportion of tertiary industry was higher than 0.1, but was not a significant driver of UR. However, the q value continued to increase from 2014 to 2022, indicating an increasing impact on the spatial distribution of UR. Simultaneously, the p value continued to decrease, suggesting that optimizing urban industrial structure has a positive impact on improving UR.

Among the indexes of urban social resilience (x4, x5, x6), the q values gradually declined from 2014 to 2022, signifying that the influence of social resilience indexes on the spatial heterogeneity of UR gradually weakened. Nevertheless, in 2022, the q values of the doctor density rose, which might be associated with the COVID-19 outbreak commencing at the end of 2019. After the COVID-19 outbreak, many doctors in China were engaged in epidemic prevention and control, facilitating each city to promptly resume normal operation. The p values of the average salary of on-the-job employees and doctor density were < 0.1, indicating that the factors x4 and x6 had a significant effect on enhancing the UR.

The q values of urban ecological resilience factors (x7, x8, x9) were small and the values of p were greater than 0.1, indicating that they had a minimal impact on the spatial heterogeneity of UR. The q values for x7 and x8 gradually increased from 2014 to 2022, although they remained low, suggesting an increasing impact on the spatial heterogeneity of UR. The q value of the sewage treatment rate was large, indicating that the control and purification of pollutants positively influenced the improvement of UR. However, this q value gradually decreased from 2014 to 2022, indicating a diminishing impact on the spatial heterogeneity of UR.

The degree of influence of driving factors on the spatial heterogeneity of UR varied significantly across different time periods. The driving factors followed the order of x4 > x1 > x3 > x6 > x5 > x9 in 2014; x3 > x4 > x1 > x6 > x5 > x9 in 2018; and x1 > x3 > x6 > x4 in 2022. The major contributing factors were the average wage of on-the-job employees and per capita GDP in 2014; per capita social consumption, average wage of on-the-job employees, and per capita GDP in 2018; and per capita GDP, per capita social consumption, and doctor density in 2022. The impact of eco-environmental factors on the spatial heterogeneity of UR was found to be small.

3.7. Detection of Interactions

The geodetector statistical tool was used to detect interactions between the driving factors of UR in Hunan Province (Table 4). In Table 4, the values on the diagonal (bold and italic) are the degrees of influence of each factor on the spatial heterogeneity of UR (q statistics); the values in the lower left corner of the diagonal are the degrees of influence of the interaction between two factors on the spatial heterogeneity of UR (q statistics); and the values in the upper right corner of the diagonal are the interaction increment of the interaction between two factors on the spatial heterogeneity of UR (white font).

For single-factor detection in 2014, the degrees of influence of q(x2), q(x7), and q(x8) with green background cells on the diagonal of Table 4 were very low, suggesting a non-significant impact on the spatial heterogeneity of UR. Other factors such as x4 had the greatest impact, followed by x1, x3, and x6 on the spatial heterogeneity of UR.

For the detection of interactions in 2014, the top three interaction detections were q(x1∩x8) = 0.9964, q(x4∩x8) = 0.9960, and q(x6∩x8) = 0.9954. The q(x2∩x7), q(x2∩x8), and q(x7∩x8) values were not significant due to their low single-factor detections (x2, x7, and x8). However, the interaction detections between these single factors and other factors was large, suggesting a significant impact on the spatiotemporal heterogeneity of UR.

In terms of the interaction increment between two factors in 2014, the positive increments indicated that factor interaction can enhance UR. The values of I_i(q)(x2∩x3), I_i(q)(x2∩x4), I_i(q)(x2∩x5), I_i(q)(x2∩x6), I_i(q)(x2∩x9), I_i(q)(x7∩x9), and I_i(q)(x8∩x9) were large, due to the nonsignificant impact of x2, x7, and x8 on UR resulting in a low degree of influence. The interaction increments between other pairs of factors were small; however, their interaction made a large contribution to UR.

For single-factor detection in 2018, the degrees of influence of the single factors (x2, x7, and x8) were comparable to those in 2014. It was found that x3 and x4 exerted the most significant impact on the spatial heterogeneity of UR, followed by x1 and x6.

In 2018, the interaction detections between two factors revealed notable combinations with high interaction detection values, specifically x2∩x3 and x2∩x4 (0.9909), x3∩x6 and x4∩x6 (0.9698), and x3∩x8 and x4∩x8 (0.9675). Conversely, the interactions q(x2∩x7), q(x2∩x8), and q(x7∩x8) were found to be not-significant, mirroring the findings from 2014.

The interaction increments between two factors in 2018 closely resembled those in 2014.

In 2022, the analysis of single-factor detections indicated that the degrees of influence of q(x2), q(x5), q(x7), q(x8), and q(x9), highlighted by green background cells on the diagonal of Table 4, were notably low, suggesting a non-significant impact on the spatial heterogeneity of UR. In contrast, factors x1, x3, and x6 had a substantial impact on the spatial heterogeneity of UR.

Regarding detection of interactions in 2022, the three highest interaction detection values were recorded as q(x1∩x5) = 0.9977, q(x1∩x2) = 0.9973, and q(x3∩x6) = 0.9938. Notably, x6 had a significant influence on UR in 2022, largely attributed to the outbreak of COVID-19. Meanwhile, the values for q(x5∩x7), q(x5∩x9), and q(x7∩x9) were relatively low. However, the interaction detections values between the single factors (x5, x7, and x9) and other factors were substantial, indicating a significant impact on the spatiotemporal heterogeneity of UR.

For the interaction increments in 2022, all factor interactions were positive, suggesting that such interactions can enhance UR. The interaction increments I_i(q)(x2∩x3), I_i(q)(x2∩x4), I_i(q)(x2∩x5), I_i(q)(x2∩x8), and I_i(q)(x2∩x9) were notably high, which was attributed to the minimal influence of x2 on UR. Conversely, the interaction increments among the other pairs of factors, expect for x2, were comparatively small; however, their interaction contributed significantly to UR.

The impact of the interaction between any two factors on the spatiotemporal heterogeneity of UR was higher than that of a single factor, suggesting that the UR of Hunan Province was influenced by multiple factors rather than a solitary one. There were two types of interaction outcomes for all factors, i.e., two-factor enhancement and nonlinear enhancement, with no linear weakening and independence between two factors observed in this study. The interaction detection value on the right side of the lower left corner of Table 4 gradually turned green over time, i.e., the detection value decreased from 2014 to 2022, indicating that the UR of ecological factors displayed a gradually weakening tendency. Conversely, the interaction detection value on the left side of the lower left corner of Table 4 increased progressively, suggesting that the influence of economic factors on UR gradually increased. This was because economic development served as the foundation for other resilience improvements and sustainable urban development.

To summarize, the interaction detections at each time node were essentially the same although the spatial heterogeneity of UR varied at different time nodes. The factors contributing to the spatiotemporal heterogeneity of UR in Hunan Province mainly arose from the disparities in economic development and social security among different cities. The coordination and joint action of various driving factors was the main reason for the spatiotemporal heterogeneity of UR.

4. Discussion

Cities are complex human-environment coupling systems and the processes by which they develop resilience are comprehensively affected the alterations of the internal and external urban environment. The spatiotemporal heterogeneity of UR is the outcome of the joint effect of various natural environment, economy, society, ecology, facilities, and internal and external interferences. Identifying the spatiotemporal heterogeneity of UR and its key determinants holds significant theoretical and practical significance for determining the UR mechanism and formulating the strategies of resilience promotion.

The conclusions of this study were largely consistent with the distribution principles of China's UR (Zhao, 2021) and Hunan’s UR (Chen, Ma, Yang, Zhou, Huang, & Chen, 2023). For example, high levels of resilience were identified in cities surrounding a provincial capital city; high-resilience cities had a positive spatial spillover effect; and low-resilience cities generally tended to develop toward high-resilience in cities over time. Hunan's UR increased from 2014 to 2022 and was high in the northeast and low in the southwest; The spatiotemporal heterogeneity of UR was mainly driven by economic and social factors.

The average value of Hunan's UR (0.2692 in 2014, 0.3021 in 2018, and 0.3422 in 2022) was lower than the average value of China's UR of 0.3604 in 2010 (Zhao, Fang, Liu, & Zhang, 2022). This indicates that Hunan's UR had not caught up with China's average resilience by 2022, and Hunan's economic, social, and ecological development remains at a low level. Direct comparison between the conclusions of Zhao et al. (2021) and the present study is difficult because Zhao et al. focused their resilience analysis on 292 cities, rather than all cities in China.

Previous studies of the driving factors of UR predominantly employed correlation and partial correlation coefficients. These methods mainly concentrated on the linear relationship between a single driving factor and the research object. In reality, the spatiotemporal heterogeneity of UR is typically a consequence of the interactions among multiple factors. The geodetector statistical tool used in this study can not only explore the impact of a single factor on UR, but also determine the contribution of factor interactions on UR.

The directionality of the spatial distribution of UR based on the standard deviation ellipse can describe the centrality, distribution, directionality, and spatial form of the spatial distribution of UR from a global and spatial perspective. The spatial range, barycenter, major axis, minor axis, and azimuth of the ellipse were used to reveal the outline and dominant direction of the spatial distribution of UR.

The distribution of the Theil index of UR in Hunan Province is associated with the imbalance of its internal economic development. As the political, economic, and cultural center of Hunan Province, Changsha has well-developed manufacturing, service, and high-tech industries due to its favorable geographical location and abundant resources, along with substantial policy support and resource prioritization. In contrast, western and southern Hunan are located in mountainous terrain, characterized by challenging transportation access, limited resource availability, and a reliance on a singular industrial structure. Therefore, achieving the balance of UR and sustainable economic and social development is the primary task confronting Hunan Province.

The following measures should be adopted to enhance the UR in Hunan Province. First, economic development constitutes the most significant factor influencing the spatial heterogeneity of UR (Davidson, Nguyen, Beilin, & Briggs, 2019). The high-quality and diversified urban economy with high levels of innovation, particularly in high-tech industry and the modern service sector, can foster new sources of economic growth. Additionally, emphasis should be placed on industrial transformation and upgrading, such as the transformation from old to new driving forces and the transformation of economic growth patterns. Concurrently, the anti-risk capability of the urban economy should be enhanced through measures including the optimal allocation of resources, industries, and skills as well as increased investment in urban infrastructure.

Second, social development serves as the foundation for the improvement of UR in Hunan Province. Urban social resilience should be enhanced through measures such as increased investment in governance resources and incomes for urban residents; the establishment of emergency management systems, public utilities, smart cities, smart transportation, and smart healthcare; and the improvement of population quality (i.e., education levels, health status, skills and competencies, workforce capabilities, social cohesion).

Finally, ecosystem resilience plays a crucial role in ensuring the smooth functioning of urban systems and creating a comfortable living environment in Hunan Province. The ability of cities to cope with various ecological risks would be enhanced through measures, such as eco-environmental protection, the development of low-carbon and green industry, resource recycling, and integrated and differentiated construction approaches.

5. Conclusions

Identifying the key determinants and their contributions to the spatiotemporal heterogeneity of UR provides the basis for ensuring the balanced development of UR, formulating urban planning, and achieving urban sustainable development. A UR evaluation model was constructed to reveal the spatiotemporal heterogeneity of UR in the GIS environment. The aggregation characteristics of high and low values of UR were revealed by applying the hotspots method. The evolution and direction of UR were spatially and dynamically determined using the standard deviation ellipse. The overall, interregional, and intra-regional differences of UR were determined by applying the Theil index. The degrees of influence of a single factor and interactions between two factors on the spatiotemporal heterogeneity of UR were determined using the geodetector tool.

From 2014 to 2022, the UR in Hunan Province generally increased; the resilience gap between cities gradually narrowed; the number of high-resilience cities gradually increased; and the resilience, hotspots, and Theil index were all high in the northeast and low in the southwest. High resilience cities had a positive spatial spillover effect on surrounding low-resilience cities. The barycenter of UR shifted from the northeast to southwest and then to southeast during the study period. The factors that contribute significantly to UR were per capita GDP, average wages of on-the-job employees, per capita social consumption, and doctor density. The impact of the interactions between any two factors on the spatiotemporal heterogeneity of UR was higher than that of a single-factor. Two-factor enhancement and nonlinear enhancement were the main forms of the factor interaction of UR in Hunan Province. Economic development and social security constitute the main material basis for the spatiotemporal heterogeneity of UR in Hunan Province.

The measurement of UR should encompass a unified theoretical framework featuring multiple elements (environment, facilities, and system), multiple levels (communities, cities, and regions), and multiple objectives (safety, adaptation, and innovation). The effectiveness of decision-making for UR should be enhanced through dynamic monitoring, simulation, and early warning systems. The mechanism of the spatiotemporal heterogeneity of UR could be determined by considering the resilience processes with long-time series and classifying urban types according to different regions, sizes, and functions. The coupling and causal relationships between urban economic, social, and ecological factors and UR should be comprehensively investigated using a Bayesian probability network, structural equation model, and the system dynamics method.