System-Reliability-Based Disaster Resilience Analysis of Infrastructure Networks and Causality-Based Importance Measure

The socioeconomic activities of modern societies, driven by rapid population growth and changes, are supported by increasingly large and complex infrastructure facilities. Under this circumstance, the escalating climate change and the interactions of components within the facilities make civil infrastructure systems more vulnerable and thus pose unprecedented risks. For this reason, research has been conducted to evaluate the impact of disasters on infrastructure systems from the viewpoint of disaster resilience. In the literature, disaster resilience is mainly defined as the ability of communities or systems to adapt to, recover from, and prepare for hazards, shocks, or stresses such as earthquakes, droughts, or violent conflicts.

Bruneau et al. [1] proposed a multi-aspect view of disaster resilience, which characterizes the target system by four properties: robustness, redundancy, resourcefulness, and rapidity. Using this concept, Chang et al. [2] developed resilience measures related to the expected loss and applied them to community performance objectives. As such, initial studies of disaster resilience and reliability of infrastructure systems were mainly performed for community-level assessment. However, the concept can be extended to various scales and types of civil infrastructure systems, e.g., lifeline networks and structural systems.

Earlier research on disaster resilience focused on calculating the system reliability of infrastructure networks, e.g., power distribution systems [3,4], pipelines [5,6], transportation networks [7–9], and bridge networks [10,11]. Researchers also attempted to identify critical risk factors and component failure combinations from a system performance degradation viewpoint [12,13]. It was found that the computational cost required for system reliability calculation grows exponentially with the number of network components. Several approaches have been developed to address the issue. Byun [14] proposed a matrix-based Bayesian network (MBN) scheme that facilitates probabilistic inference on complex systems for given evidence through a novel data structure storing and handling information regarding component reliability. Approximate methods based on alternative system representation or calculations were also developed [15–18].

Further research efforts incorporated the two directions discussed above into the pre- and postdisaster stages of network resilience analysis. Jönsson et al. [19] and Johansson et al. [20] performed reliability and vulnerability analyses for power networks to demonstrate the critical need for system-level analyses in disaster resilience assessment. The vulnerability analysis aimed to quantify the negative system-level consequences for given characteristics of disasters or the combinations of component failures. Kilanitis and Sextos [21] performed the fragility analysis of each road network component based on seismic fault information and investigated changes in network functionality and loss factors caused by traffic delay using the samples of the component states. Considering the impact of network component capacities, Zhang and Alipour [22] optimized the risk mitigation strategies of components in the predisaster stage to minimize direct and indirect costs of the whole network in the postdisaster stage.

Beyond evaluating the disaster resilience of networks, various attempts are increasingly made for resilience-informed decision-making regarding network design and maintenance strategies. To keep up with such demands, we should be able to prioritize network components in terms of their contributions to system resilience. For this purpose, various importance measures (IMs) were proposed to quantify the influences and contributions of the components [23]. Probabilistic interpretations of widely used IMs were categorized in terms of the information used for calculating IMs. Beyond the traditional definitions, IMs were further developed to reflect various perspectives and conditions, e.g., topological viewpoint and reliability [24], and multiple states of components [25,26]. A well-defined IM can provide insights for component-level decision-making in many practical problems, such as finding the priority of reinforcement or inspection [27,28].

To serve as an informative measure from a disaster-resilience perspective, IM should be able to describe the impact of the component states on the system-level performance considering intercomponent correlations [29]. However, it is noted that the existing IMs in the literature may not successfully exclude the effects of initial disruption scenarios making insignificant contributions to the system performance degradation. Moreover, the evaluations of components' relative importance should consider both the reliability of individual components and the physical impact of component failures on the system's performance.

Therefore, we apply the system-reliability-based disaster resilience analysis framework [30] to infrastructure networks for the first time to facilitate resilience-informed decision-making processes of the systems. Then, we propose a new IM based on the resilience analysis results and a causal model of the system-level failure. The proposed causality-based importance measure handles network failure through major disruption scenarios. Without losing applicability to general network reliability problems, this paper focuses on those in which the system performance is defined by the connectivity between major components, which is a top priority to secure evacuation routes or production procurement paths in the postdisaster stage.

The rest of the paper is organized as follows: Section 2 briefly reviews the system-reliability-based disaster resilience analysis framework [30], widely known IMs in the literature, and theories of causality effect evaluation by a causal diagram. Section 3 provides a hypothetical network to demonstrate the general applicability of the framework to infrastructure network scale and IM calculations. Section 4 proposes new causality-based IMs based on the disaster resilience analysis framework and a causal model. The numerical examples in Sec. 5 further investigate the effects of the correlations between component failures and test the applicability and merits of the proposed causality-based IM. The paper concludes with a summary and concluding remarks in Sec. 6.

Lim et al. [30] pointed out the limitations of existing disaster resilience analysis methodologies: the definition of system performance measure in a restoration curve is often subjective or elusive; the interaction between the system and its components is not incorporated well into the analysis; and the methodologies often lack a systematic procedure to utilize the analysis results as a basis for the decision-making process to secure the disaster resilience of the system. Hence, a system-reliability-based disaster resilience analysis (S-DRA) framework was proposed to overcome the limitations.

The S-DRA framework features three criteria of disaster resilience — reliability, redundancy, and recoverability. Reliability is the component-level capability to avoid or minimize initial component disruptions during disastrous events. Redundancy is the capability of a system to minimize system-level performance degradation due to cascading failure initiated by component-level disruptions. Recoverability is the ability of engineers and society to take appropriate actions on components and systems to restore the degraded system's functioning quickly and wholly. The three criteria are illustrated in Fig. 1.

The three criteria of disaster resilience should be defined and evaluated considering the scale-specific characteristics because system modeling, analysis of disaster-induced damage, and recovery strategy all vary depending on the scale of the analysis. In this context, Lim et al. [30] laid out a “3 × 3 resilience matrix” whose rows and columns, respectively, describe the three analysis scales and the three criteria, as shown in Table 1. Each element in the resilience matrix provides detailed descriptions of the resilience criterion for the given analysis scale.

While Lim et al. [30] focused on the individual structure scale, i.e., the first row of the matrix, further developing the framework and proposed detailed analysis methods based on single structure examples. This paper rather focuses on the second row of the resilience matrix to facilitate disaster resilience evaluations and decision-making for infrastructure networks such as transport networks, power grids, hydraulic systems, and gas distribution systems. It is noted that such networks are often modeled as graphs consisting of nodes representing bridges and stations and links such as transmission lines and pipelines. The nodes and links in a graph model are treated as node- and line-type components of the target network system.

The reliability evaluation at the network scale should properly model and consider the variability of the hazard intensity, the distribution of network components over a large area, and the spatial correlation of hazard demands on components. In evaluating the redundancy of a network system, the effects of sequential component failures triggered by initial component disruptions on the eventual network-level performance degradation should be considered. For recoverability evaluations, we should be able to assess the overall ability of engineers and society to recover the damaged network in the postdisaster stage and secure essential resources and response manuals in the predisaster stage [31]. Following the definition of the criterion, diverse socioeconomic factors describing the network environment should be studied carefully. Therefore, this paper focuses on infrastructure networks' reliability and redundancy criteria (the shaded areas in Table 1) while leaving recoverability as a future research topic.

To facilitate S-DRA exploiting the reliability index $(β)$ and redundancy index $(π),$ a reliability-redundancy (⁠$β$-$π$⁠) diagram was proposed. The axes of the two-dimensional scatter plot represent the redundancy and reliability indices, respectively. A point with the coordinates $(πi,βi)$ in a $β$-$π$ diagram represents the indices computed by Eqs. (1) and (2) for an initial disruption $Fi$ under the given hazard $H.$ The $β$-$π$ diagram serves as a tool to provide a holistic description of the disaster resilience of the system both at the component level and the system level and to support relevant decision-making processes.

where $Φ(⋅)$ is the standard normal CDF and $Pdm/λH=Hdm$ was termed “per hazard de minimis risk” of the hazard $H.$ The boundary of the domain $DλH$ serves as the disaster resilience limit-state surface in the $β‐π$ diagram, as shown by the contours of the $β‐π$ diagrams in Ref. [30] and this paper. The disaster resilience analysis procedure employing the $β‐π$ diagram was termed a “reliability-redundancy $(β‐π)$ analysis.”

In risk management of complex systems, quantifying the components' relative contributions to the likelihood of the system's failure is often helpful. To this end, importance measures (IMs) have been developed to guide risk-informed decision-making regarding various systems. We briefly review the following IMs, widely used in system reliability analyses [33,34]. These IMs will be demonstrated through the hypothetical network example (Sec. 3) and compared with the causality-based IM proposed in this paper (Sec. 4).

where $Fsys−e$ and $Fsys+e$ respectively denote the failure events of the target system with component $e$ removed (or failed) and replaced by a perfectly reliable component.

The IMs in Eqs. (5)–(10) provide larger values for components contributing more to the system failure probability. By definition, RAW and RRW are greater than or equal to 1, while the other IMs range between 0 and 1. In calculating RAW, RRW, and BP, the state of the component $e$ is fixed regardless of their likelihood of failure.

A causal model aims to provide mathematical descriptions of the causal relationship between random variables. In particular, recent research has often employed a causal diagram, i.e., a directed acyclic graph consisting of the nodes representing random variables and the arrows describing the causal influence between those variables [36,37]. A causal diagram well constructed using relevant variables can play an essential role in answering causality-related queries such as “How much would the disruptions of certain components affect the network's performance?” based on observable data. (Controlled experiments for statistical causality identification are impossible for already-built civil infrastructure systems. Therefore, this study utilizes computational simulation data obtained for each component disruption scenario. Then we process the data using causal diagrams to consider causal effects in quantifying the relative importance of components.)

It is noted that direct inference based on observable data may lead to spurious misunderstanding of causality due to correlations. For example, let us consider a causal diagram in Fig. 2(a) that describes the relationship between binary random variables. It is noted that $S$ and $Z$ may represent a set of random variables. Suppose the conditional probability $P(Y=yX=x)$ is computed using Bayes' rule to answer the query, “How much does $X$ affect $Y$?” Although there are no directed paths from $X$ to $Y,$ a probabilistic inference that may give a result $P(Y=yX=x)≠P(Y=yX≠x).$ Such statistical dependence between $X$ and $Y$ is due to the common source effects from $S,$ not the causal relationship between $X$ and $Y.$

where $do(X=x)$ is an operator representing the intervention shown in Fig. 2(b). The causal effect can be quantified by the change in the conditional probability caused by the opposite outcome, $X≠x.$ For example, the difference $P(Y=ydo(X=x))−P(Y=ydo(X≠x))$ or the ratio $P(Y=ydo(X=x))/P(Y=ydo(X≠x))$ are often adopted depending on the analysis purpose [39].

Various formulaic test methods have been developed to identify which variables should be observed and controlled to investigate the causal relationship of interest. Among those methods, this study adopts a simple graphical test called a “back-door criterion.” In general, a variable set $A$ is referred to as meeting the back-door criterion for a variable pair $(X,Y)$ in a directed acyclic graph $G$ if they satisfy the following conditions: (1) no elements in $A$ are descendants of $X$⁠, and (2) $A$ blocks every path between $X$ and $Y$ that contains an arrow into $X.$ The name of the criterion, back door, is attributed to the second condition in which the path with arrows heading to $X$ is interpreted as entering $X$ through the back door. The phrase “$A$ blocks every path” in the second condition means that $X$ and $Y$ become independent by the observation of the variables in $A.$ (For example, for the causal diagram in Fig. 2(a), $S$ and $Z$ satisfy the back-door criterion for a variable pair $(X,Y)$ because they block all back doors from $X$ to $Y,$ along the path $X–S–Z–Y.$⁠)

because the state of $Y$ depends only on the state of $Z.$ The causal query $P(Y=y|do(X≠x))$ can be computed in a similar way using the back-door criterion in Eq. (12), which leads to the same value as $P(Y=y|do(X=x)).$ This consistency on $Y$ indicates that there exists no direct causal effect of $X$ on $Y,$ as the causal diagram in Fig. 2(a) shows.

This section presents a new approach from the basic modeling of infrastructure networks to the creation of initial disruption scenarios and the subsequent realization of a resilience index. Let us consider an infrastructure network consisting of $n$ components with binary states – failure or nonfailure. All combinations made by $k$ components are represented by the component index sets$Iik, i=1,…,(nk)$⁠. For instance, when $n=5$ and $k=2$⁠, the index sets are $I12={1,2},$$I22={1,3},$…, $I102={4,5}.$ For each combination, a mutually exclusive and collectively exhaustive (MECE) scenario set$Fik$ is introduced to enumerate all joint states. For example, $F12={E1E2, E1E¯2, E¯1E2,E¯1E¯2}.$ In particular, the corresponding disruption scenario is denoted by $Fik.$ For example, $F12=E1E2.$

where $Xe$ is a Bernoulli random variable indicating the failure by the value 1 and nonfailure by the value 0. The reliability of $n$ components can be represented by an $n$-variate standard normal distribution $Z∼Nn(0,Σ)$ with a mean vector of $0∈Rn$ and a covariance matrix $Σ,$ which is equal to the correlation coefficient matrix $R∈Rn×n.$ The correlation coefficient between different $Ze$'s can be computed based on the marginal and joint failure probabilities [17].

The probability $P(⋃l(Cl∩e∈IikEe)H)$ can be computed by system reliability analysis using the matrix-based system reliability method [29], sequential compounding method [16], or Monte Carlo simulation.

Figure 3(a) shows the configuration of a hypothetical infrastructure network consisting of nine node-type components near the seismic source. Suppose the reliability indices of the components $βe, e=1,…, 9$ are evaluated as {1.600, 1.718, 2.014, 2.100, 2.181, 2.403, 2.600, 2.662, 2.836}. The covariance matrix of $Z∼N9(0,Σ)$ is determined based on the assumption that the correlation coefficient between $Ze$'s of two components with distance $Δ$ (km) is $exp(−Δ−0.5).$ The system failure $Fsys$ is defined as the disconnection of all terminal nodes from the source nodes, described by the six cut sets as $Fsys=⋃l=16Cl=E2E3∪E2E6∪E3E5∪E5E6∪E7E9∪E8.$

The redundancy index is calculated using Eq. (16) as $π12=−Φ−1(0.1849)=0.8969.$

The β-π diagram in Fig. 3(b) is obtained by repeating this process for all scenarios in $Fik (k=1, 2, 3).$ The markers in the diagram are distinguished with different shapes and colors according to $k,$ the number of components in the disruption scenario. Especially for $k=1,$ nine components are labeled near the marker with the corresponding color. The contours in the figure represent the boundaries of the disaster resilience constraints $DλH$ in Eq. (4) for the four annual mean occurrence rates $λH$ and $Pdm=10−7/year.$

The redundancy indices of the initial disruption scenarios that include one or more minimum cut sets are $−∞$ because the conditional probability $P(Fsys|Fik,H)$ is 1. For convenience, the index pairs of such scenarios are located in $π=−3$ in the diagram. This visualization is acceptable because the vertical location of the marker is almost consistent for $π≤−3,$ and thus the reliability index determines whether the system satisfies the resilience constraint.

The $β‐π$ diagram helps us to identify critical initial disruption scenarios leading to socially unacceptable risk. For example, when $λH=1/4,800year,$ among the nine initial disruption scenarios with $k=1,$ only $E7$ and $E9$ satisfy the disaster resilience constraint $DλH$ in Eq. (4). It is observed that the markers tend to move in the upper left direction as $k$ increases, i.e., higher reliability index and lower redundancy index. This is a natural phenomenon because the spatial demand on the components and correlation coefficients between component states are driven by the seismic hazard with the specific mean occurrence rate. In detail, stemming from the given hazard, the likelihoods of scenarios involving more components are relatively low, but once such disruption occurs, the network's capability to avoid system-level failures drops significantly. Incorporating the information on recoverability into the diagram by the colors or sizes of the markers [30], the $β‐π$ diagram can support decision-making processes to manage the network's disaster resilience.

Next, the six IMs reviewed in Sec. 2.3 are calculated for the hypothetical network, as shown in Fig. 4. Six IMs quantify the relative importance of components differently based on the reliability of individual components and their topological importance in the network. The existing IMs commonly find components 1, 4, 7, and 9 relatively less critical but provide different comparisons for the others. It is also noted that these IMs do not check the contributions of the components to inducing socially unacceptable risks while considering the hazard occurrence rate. Besides, the CP and ICP of components 1 and 4 are nonzero even though they do not contribute to the system performance from the network topology viewpoint (because the state of component five governs). Further in-depth analysis of these IMs and comparison with the causality-based IMs will be discussed in Sec. 4.2.

To reflect the causal effects of component failures on the system failure properly, first, we describe the causal relationship between various factors, e.g., hazards and component states, using a causal diagram. Then, from the scenario-wise S-DRA, we allocate the results to individual components considering the causal contribution calculated through the causal diagram and propose new IM. Finally, we discuss and compare the resilience philosophy of the novel causality-based IM with the existing correlation-based IMs and do-free IMs set as the comparison group.

Figure 5(a) shows the proposed structure of a causal diagram $G$ describing the Bernoulli random variable $Ynet$ representing the network connectivity loss induced by external hazards denoted by $S1, S2,$ …, $Sm.$ The node $U$ represents a set of unobservable factors affecting $n$ component failure variables $X1, X2,$ …, $Xn. U$ remains an implicit variable acknowledging the imperfect knowledge of the component failures and their correlations.

The node $U$ and the $m$ hazard sources are common causes of the component failure $Xe,$ which may eventually lead to network connectivity loss $Ynet.$ These common causes are the confounders inducing associations between component failures and those between component failures and a network connectivity loss, e.g., $X1→U→X3,$ and $X1→S1→X2→Ynet.$ These associations are considered spurious because the observation of the failure of a component does not necessarily change the failure mechanism of other components from a causal viewpoint.

The causal effect on the left-hand side can be estimated by computing all the probability terms using multivariate normal CDFs.

In $β‐π$ analysis, the disaster resilience of a system is evaluated for each disruption scenario, which is described as a single component failure or the joint failure of multiple components. The causal contributions of individual components within a multi-component-failure scenario are identified through the following process. First, to deal with the effect of change in the state of component $e∈Iik,$ the remaining $(k−1)$ components are denoted by the sorted index set$Ii,ek=Iik⃥{e}.$ As $Fik$ was enumerated for $Iik,$ an MECE scenario set $Fi,ek$ can be obtained for $Ii,ek.$ For example, let us consider the component index set $I13={1, 2, 3}.$ When we investigate the importance of component $e=2, I1, 23=I13⃥{2}={1, 3}$ and $F1, 23={E1E3, E1E¯3, E¯1E3, E¯1E¯3},$ respectively. For $k=1, Fi, ek$ is defined as a universal set $U$ for the given sample space.

Table 2 shows five cases of the possible values of the probability terms, along with the corresponding causal contribution of component $e, rek(F)$ computed by Eq. (20). In Case 1, the value of the conditional probability $P(FsysF,H)$ is set to zero, i.e., $F$ includes a link set. On the other hand, in Case 5, the conditional probability is set to 1, i.e., $F$ includes a cut set. In these cases, the system's state is determined regardless of the state of component $e.$ Therefore, it is found reasonable that $rek(F)$ takes zero in these cases. Case 2 indicates that the nonfailure of component $e$ guarantees the system-level nonfailure, which makes $rek(F)=1$ from Eq. (20) reasonable. Case 3 is a general case in which the contribution of the changed state of component $e$ determines $rek(F).$ In Case 4, the change in the component's state considered by the do-operator leads to no changes in the conditional probability, and naturally, Eq. (20) gives $rek(F)=0.$

Using Eqs. (22)–(24), we calculated NCIMs of the hypothetical network in Sec. 3.2 for the seismic hazard with $λH=1/4,800year$ and $1/1,000year,$ as shown in Figs. 6(a) and 6(b), respectively. It is noted that component 8 shows the most significant decreases in its NCIMs as the number of initially disrupted components (⁠$k$⁠) increases because the component contributes to the system-level failure through its topological importance rather than through joint failures with other components. Similar to the existing IMs in Fig. 4, the NCIMs of components 1, 4, 7, and 9 are insignificant.

An essential distinction of the proposed NCIM from the existing IMs is that it can evaluate the component importance for the given occurrence rate of the hazard by the $β‐π$ analysis. For example, when $k=3,$ the NCIMs of components 5, 6, and 8 increase significantly as $λH$ increases. It is also noteworthy that NCIM can consider the causal effects. For instance, the NCIMs of components 1 and 4, making no topological impacts on the network connectivity according to Fig. 3(a), are zeros. This is because the do-operators filter the spurious effect of the correlations in Eq. (20). By contrast, the CP and ICP, defined as conditional probabilities, show non-zeros values for components 1 and 4 as shown in Figs. 4(a) and 4(b). Therefore, CP and ICP cannot exclude spurious association from component failures to system-level failure, and ICP reflects only the redundancy aspects of the $β‐π$ analysis.

Components 7 and 9 are located far from the seismic hazard to have high reliability; thus, the corresponding FVs and RRWs are close to the minimum value. However, the definitions of FV and RRW fundamentally make it impossible to distinguish the relative importance of the components belonging to the same cut set. In summary, the existing IMs defined as a single ratio or difference between two probabilities have limitations in reflecting the reliability and redundancy perspectives in contrast to the proposed NCIM.

For the hypothetical network, by following the process of Eqs. (22)–(24) using the correlation-based contribution in Eq. (30) instead of Eq. (20), the normalized do-free importance measure (NDIM) is calculated as shown in Figure 7 for comparison.

In addition, NDIMs assign a higher importance weight to component 8 as $k$ increases, whereas the weight decreases when NCIM is used. It implies that the do-operator helps NCIM avoid placing extra weights on components contributing to the network redundancy for higher $k$⁠.

A numerical example is created using the topology of a bridge network in Sioux Falls, to test the applicability of the S-DRA framework and the causality-based IM to general infrastructure networks. The reliability and redundancy indices are computed for bridge failure scenarios based on a seismic hazard model and the fragility analyses of individual bridges. The impacts of the correlation coefficient between the component state on the $β‐π$ analysis process and NCIM calculation are investigated in detail. The resilience characteristics of the network under multiple seismic hazards are also examined.

Figure 8(a) shows the road map of the Sioux Falls area, where the highways and the lower-level roads are located along the outskirts and the region's interior, respectively. Figure 8(b) shows the bridge network model consisting of 19 components with two types of bridges, i.e., multispan simply supported concrete girder (MSC) for the highway bridges and single-span concrete girder (SSC) for the bridges on downtown roads. The network-level performance goal is the posthazard connectivity between the two three-node sets, representing the source and terminal of the traffic.

where $Dje$ denotes the seismic demand intensity, $F(⋅)$ is a ground-motion model (GMM), $Mj$ is the earthquake magnitude, $Rje$ is the distance between the epicenter and the structure, and $ηj$ and $εje$ are inter- and intra-event uncertainties (residuals), respectively [40]. These residuals are assumed to be statistically independent and follow the zero-mean normal distributions. Therefore, if $ση2$ and $σε2$ denote the variances of the residuals, the total variance of $lnDje,$ i.e., $σT2$ is $ση2+σε2.$

where $ρe1e2(Δ)$ is a spatial correlation model, given as a function of the distance $Δ$ between two components $e1$ and $e2$ [8]. The spatial correlation starts from 1 at $Δ = 0,$ decreases monotonically, and converges to 0 as $Δ$ goes to infinity.

where $e1, c1, c2, c3, Mref, Rref,$ and $h$ are the model parameters whose values are summarized in Table 3, along with the standard deviations of the two residuals. The GMM in Eq. (36) is for $VS30=$ 760 m/s where $VS30$ denotes the time-averaged shear-wave velocity over the top 30 m near the unspecified fault. The spatial correlation is assumed to be $ρe1e2(Δ)=exp(−0.27Δ0.4)$ [17]. In addition, the structural parameters of the bridges with MSC and SSC for extensive $(EX)$ damage state are obtained as $(αMSCEX,βMSCEX)=(0.83, 0.65)$ and $(αSSCEX, βSSCEX)=(2.62, 0.90),$ respectively from the prior structural analysis [43]. For a 7.5-magnitude western earthquake $EQwest$ whose epicenter is located at $(−4,6)$ in Fig. 8(b), the failure probability of the bridges with their types and the correlation coefficient between the bridge failures are calculated using Eqs. (34)–(36). Figure 9(a) shows the failure probabilities of the bridges. Figure 9(b) shows the positive correlation coefficients between failures of the different bridges, estimated with values between 0.178 and 0.325.