Modeling the Evolution of Major Storm-Disaster-Induced Accidents in the Offshore Oil and Gas Industry

Storm disasters are the most common cause of accidents in offshore oil and gas industries. To prevent accidents resulting from storms, it is vital to analyze accident propagation and to learn about accident mechanism from previous accidents. In this paper, a novel risk analysis framework is proposed for systematically identifying and analyzing the evolution of accident causes. First, accident causal factors are identified and coded based on grounded theory (GT). Then, decision making trial and evaluation laboratory (DEMATEL) is integrated with interpretative structural modeling (ISM) to establish accident evolution hierarchy. Finally, complex networks (CN) are developed to analyze the evolution process of accidents. Compared to reported works, the contribution is threefold: (1) the demand for expert knowledge and personnel subjective influence are reduced through the data induction of accident cases; (2) the method of establishing influence matrix and interaction matrix is improved according to the accident frequency analysis; (3) a hybrid algorithm that can calculate multiple shortest paths of accident evolution under the same node pair is proposed. This method provides a new idea for step-by-step assessment of the accident evolution process, which weakens the subjectivity of traditional methods and achieves quantitative assessment of the importance of accident evolution nodes. The proposed method is demonstrated and validated by a case study of major offshore oil and gas industry accidents caused by storm disasters. Results show that there are five key nodes and five critical paths in the process of accident evolution. Through targeted prevention and control of these nodes and paths, the average shortest path length of the accident evolution network is increased by 35.19%, and the maximum global efficiency decreases by 20.12%. This indicates that the proposed method has broad applicability and can effectively reduce operational risk, so that it can guide actual offshore oil and gas operations during storm disasters.


Introduction
Offshore oil and gas production faces a number of harsh environments. In recent years, the number of weather-related natural disasters has grown sharply as have the costs associated with related losses [1,2]. Storm disasters, such as typhoons, hurricanes, tornadoes, and storm surges, significantly impact offshore oil and gas production, and can cause operational shutdown, facility damage, casualties, ecological crises, and even serious social unrest [3][4][5]. During Hurricane Harvey (2017), about 25% of the oil and gas production in the Gulf of Mexico was shut down, and 105 production platforms (15% of the total) were closed, causing a fluctuation in oil and gas prices [6]. To avoid the repeated occurrence of accidents, lessons should be drawn from previous accidents [7][8][9][10][11]. Therefore, it is necessary to analyze the causal factors of major accidents in the offshore oil and gas industry that are caused by storm disasters and explore the accident development process. 2 of 27 To analyze the accident propagation path, it is necessary to carry out an accident investigation and identify causal factors or risk sources. At present, the most commonly used methods for accident analysis in offshore oil and gas operations include RCA, FTA, ETA, barrier analysis, and AcciMap [12]. RCA can identify the root causes of accidents through continuous iteration and improvement [13][14][15]. To ensure the accuracy and integrity of results, RCA requires a certain amount of judgement ability and a relatively exhaustive list of causal factors in advance. FTA is a method for graphically representing the possible causes of top events (i.e., accidents) and their causal relationship through top-down deductive analysis; however, this method cannot consider nonlinear interactions among causal factors [16,17]. ETA is primarily a proactive risk analysis method used for identifying the possible sequence of events after the initial event, but the choice of development paths largely depends on the level of individual knowledge and personal experience [18][19][20]. Through barrier analysis, the risks associated with cascading accidents can be identified to confirm the placement and performance of barriers. Combined with other methods (such as ETA), preventive measures for cascading accidents can be identified for offshore oil and gas operations [21,22]. AcciMap is not only a pure accident investigation tool, but can also capture socio-technical factors and illustrate the interactions among these factors causing corresponding events. With AcciMap, human and organizational factors in offshore oil and gas accidents can be highlighted [23][24][25]. All these methods have advantages but are highly subjective and deductive, thus requiring a strong level of expert knowledge and experience. For the same accident case, researchers with different experience and knowledge levels may obtain completely different results. This is very unfavorable in ensuring the accuracy of identifying the cause factors of accidents.
Complex systems often fail in complex ways, and accidents are usually the result of the interaction of multiple factors. Many methods have been proposed to analyze interactions among factors in accident propagation evaluation, including AHP, ANP, SEM, DEMATEL, and ISM [26,27]. AHP provides a convenient method for multi-objective and multi-criteria decision problems, which derives priority using a nominal scale through paired comparison of elements at the same level. As an extended form of AHP, ANP considers the dependence of elements among different hierarchies; however, the hierarchical structure of both methods is based on subjective expert judgments [28][29][30]. SEM measures the influence of each cause upon the effect, and the regression coefficient can be calculated to express the causal relationships among latent variables (i.e., factors) based on a large number of questionnaires [31][32][33]. By building a structural model containing causal relationships among complex factors, DEMATEL can analyze the interdependences among factors but is limited due to little information [34,35]. ISM is also an effective tool for identifying causal relationships among complex factors, and can help to clarify the hierarchy and priority of factors. Causal relationships are easy to grasp by ISM, but this method requires notable matrix computation resources [36,37]. All the above-mentioned methods have their own application scenarios and functional characteristics. A large number of questionnaires and computing resources are needed to obtain the quantitative interactions among different factors. At the same time, limited information can only judge whether factors are independent. Therefore, a single method cannot achieve efficient and accurate analysis of the interactions among factors in accident propagation evaluation.
To analyze the evolution of accidents, it is necessary to evaluate the propagation path of causative factors. BN has been widely applied for quantitative risk assessment, as forward predictive analysis and backward reasoning diagnosis can be carried out in BN [38,39]. The primary advantage of BN is its probability updating and continuous learning ability. In the field of the offshore oil and gas industry, BN is mainly used to evaluate the dynamic risk evolution of major accidents including blowouts and explosions [40,41], as well as the failure probability and reliability of key equipment and systems; examples are blowout preventer [42], submarine oil and gas pipeline [43,44], drilling riser [45], crude oil separation system [46], and managed pressure drilling system [47,48]. However, it is not easy to obtain the prior probability for BN. Besides, BN focuses on the assessment of the overall risk level of the network, and cannot deal well with the problems of local risk and single path risk. Although BN can determine the importance of nodes through sensitivity analysis, it is hard to study the importance of different paths using the method of path attack. In addition, if the number of causal factors is large and the interaction relationships are complex, there may be a state explosion problem [49].
The limitations of these methods including strong subjectivity, unclear clarification of interactions among factors, and difficulty in assessing path risks. To address issues in the above-mentioned methods, this paper proposes a novel risk analysis framework for systematic identification and evaluation accident propagation. We start with the identification of the accident causes based on GT. Then DEMATEL is integrated with ISM to analyze interactions among factors. Finally, HADY is incorporated into CN to describe the evolution process of accidents and calculate shortest paths. Compared to reported works, the detailed contributions of the proposed method are summarized as follows: 1. GT is developed to objectively identify accident causal factors from original data, which can reduce personnel subjective influences imposed by their knowledge and skill level. 2. The use of DEMATEL-ISM can compensate for the shortcomings of insufficient information and high computing costs. In addition, the way to establish the direct influence matrix is improved according to the interaction matrix obtained from the objective analysis in the previous step. 3. A hybrid algorithm that can calculate multiple shortest paths of the accident evolution under the same node pair is proposed. Combined with the causal factors and interaction relationships obtained by GT and DETAMEL-ISM, CN can quantitatively describe the evolution process of accidents.

Methodology
To analyze the evolution process of major accidents in the offshore oil and gas industry under storm disasters, a method for the systematic identification and analysis of the evolution process of accident causes is proposed and applied in a case study, which includes four main steps: 1. Identification of accident causal factors. 2. Analysis of the hierarchical structure of accident causes. 3. Research on the evolution process of accident causes. 4. Study of the application in a specific case Firstly, accident causal factors are identified and coded based on GT, and interaction relationships among causal factors are analyzed according to the occurrence frequency of causal factors in the same case. Then, using the DEMATEL-ISM method, the influencing relationships among causal factors are analyzed. When the hierarchy of causal factors is divided, a hierarchy of accident evolution is established. Subsequently, the method of CN is used to determine the evolution model of accident causes, analyze the characteristics of accident evolution, and calculate the shortest paths of accident evolution. Finally, the systematic identification and analysis of the evolution process of accident causes is applied in a case study. The specific process is depicted in Figure 1.

Identification of Accident Causal Factors
As an inductive qualitative research method rooted in original data, GT is an effective tool to identify accident causal factors. So far, there have been three schools: classic GT [50], programmatic GT [51] and constructivist GT [52]. As an objective induction method, the obtained analysis results are more objective and reliable as they are based on accidentrelated data [53]. Considering the standardization and simplicity of programmed language and in line with the method, this paper is based on programmatic GT. The coding process of programmatic GT includes three main steps: open coding, axial coding, and selective coding [54]. The method of GT emphasizes constant comparison and abstraction to form a reliable theory; thus, theoretical saturation tests should be performed after coding.

Identification of Accident Causal Factors
As an inductive qualitative research method rooted in original data, GT is an effective tool to identify accident causal factors. So far, there have been three schools: classic GT [50], programmatic GT [51] and constructivist GT [52]. As an objective induction method, the obtained analysis results are more objective and reliable as they are based on accidentrelated data [53]. Considering the standardization and simplicity of programmed language and in line with the method, this paper is based on programmatic GT. The coding process of programmatic GT includes three main steps: open coding, axial coding, and selective coding [54]. The method of GT emphasizes constant comparison and abstraction to form a reliable theory; thus, theoretical saturation tests should be performed after coding.

Analysis of the Hierarchical Structure of Accident Causes
To analyze complex interaction relationships among accident causes, a suitable and efficient analysis method is needed. As important tools for system analysis and decision making, DEMATEL [55,56] and ISM [57] have been widely used in various fields. Nevertheless, DEMATEL fails to consider the influence of factors on itself, while ISM fails to consider the strength of influence relationships among factors. The integrated DEMATEL-ISM method compensates for the deficiencies of both tools, and the integrated method also reduces the difficulty for effectively calculating the reachability matrix in ISM. By setting the threshold, DEMATEL-ISM eliminates the weak influencing relationship in the system and simplifies the structure of the system.
In this paper, the traditional method of DEMATEL-ISM is improved, by incorporating the following points: 1. When constructing the index system of causal factors, the categories obtained from GT are absorbed, and are put into the set Y (yi&yj ∈ Y, i = 1, 2, …, n; j = 1, 2, …, n) of causal factors, to reduce the influence of subjectivity and avoid the inconsistency of the index scope. 2. When establishing the matrix F, the matrix R is referred based on expert scoring, to obtain more objective analysis results.

Analysis of the Hierarchical Structure of Accident Causes
To analyze complex interaction relationships among accident causes, a suitable and efficient analysis method is needed. As important tools for system analysis and decision making, DEMATEL [55,56] and ISM [57] have been widely used in various fields. Nevertheless, DEMATEL fails to consider the influence of factors on itself, while ISM fails to consider the strength of influence relationships among factors. The integrated DEMATEL-ISM method compensates for the deficiencies of both tools, and the integrated method also reduces the difficulty for effectively calculating the reachability matrix in ISM. By setting the threshold, DEMATEL-ISM eliminates the weak influencing relationship in the system and simplifies the structure of the system.
In this paper, the traditional method of DEMATEL-ISM is improved, by incorporating the following points: 1. When constructing the index system of causal factors, the categories obtained from GT are absorbed, and are put into the set Y (y i &y j ∈ Y, i = 1, 2, . . . , n; j = 1, 2, . . . , n) of causal factors, to reduce the influence of subjectivity and avoid the inconsistency of the index scope. 2. When establishing the matrix F, the matrix R is referred based on expert scoring, to obtain more objective analysis results. 3. In the evolutionary hierarchy of accidents obtained from ISM analysis, the category frequency and the strength of influence relationships among factors obtained from GT and DEMATEL are presented simultaneously.
The specific algorithm steps are as follows: 1. Establish the matrix F = [f ij ] n×n , where f ij represents the direct influence of factor y i on factor y j (y i &y j ∈ Y, i = 1, 2, . . . , n; j = 1, 2, . . . , n. Y is a set of system factors.). {0, 1, 2, 3} respectively represent {no influence, weak influence, medium influence, strong influence} in the matrix F. 2. Normalize the matrix F and obtain the matrix C = [c ij ] n×n , where c ij ∈ [0, 1]: 3.
Develop the matrix T = [t ij ] n×n which is established to couple the direct and indirect influence relationships among factors:

4.
Calculate the b i , d i , c i and a i of each factor:

5.
Develop the matrix H = [h ij ] n×n mainly for the purpose of considering the influence of factors on itself: 6.
Develop the matrix K = [k ij ] n×n , which is established to simplify the weak influence relationship between factors and highlight the hierarchy of the system by selecting an appropriate threshold:

7.
Calculate the Q i and P i : 8. Establish the condition for hierarchical partitioning: If the above condition is satisfied, it is proved that the corresponding factors in P i can all be found in Q i . Therefore, these factors are at a higher level. At the same time, the corresponding row i and column j are deleted from the matrix K, and then step 7 and step 8 are repeated until all factors are deleted. The hierarchical structure of all factors is determined according to the order in which the factors are deleted.

Research on the Evolution Process of Accident Causes
Based on graph theory and statistical mechanics, CN is an important tool for analyzing the structural characteristics and development process of complex systems. CN can describe the process of risk evolution and effectively evaluate dependent relationships among different factors. Many complex systems or processes in nature can be described by CN [58]. CN has been widely used in a variety of fields, such as transportation networks [59,60], power systems [61,62], natural disasters [63], and offshore oil and gas production [49,58].
According to the organization form of network nodes, the four basic models of regular networks, random networks, small-world networks, and scale-free networks can be classified. The small-world network is characterized by a large clustering coefficient and small average path length. The scale-free network reflects node growth and preference dependence based on small-world network characteristics, which are closer to the characteristics of network models in the real world [64,65].

Modeling of CN
The CN graph is a data structure composed of nodes, edges, and weights. The CN graph is usually represented mathematically by G = (V, E, W). The node in V is denoted as v i by its order i. The edge linking v i with v j in E is denoted as e ij . If e ij and e ji represent the same edge, the graph is an undirected graph; otherwise, it is a directed graph. Each element in W expresses the weight of e ij , which is denoted as w ij . If each w ij is equal to 1, the graph is an unweighted graph; otherwise, it is a weighted graph. In this paper, the graph is a weighted network graph. The principal assumptions for the modelling can be summarized as follows: different events in the process of accident evolution can be represented by nodes, and the associations among events can be represented by weighted edges. Causal factors (nodes) are obtained by GT, and influence relationships (weighted edges) among factors obtained by DEMATEL-ISM. The accident evolution model is constructed by CN.

Characterization of CN
Every complex network has its unique topological structure and connectivity, and the transmission process of internal information cannot be obtained only by means of observation. Therefore, we need the specific measurement method to evaluate the characteristics of the complex network, so that we can have a deeper understanding of it. Through the discrimination and analysis of network characteristics, the critical nodes, critical paths, propagation modes, and development trends of accident evolution can be obtained. These findings can provide references for accident prevention, control and emergency. The principal models and the presentation of the main measurements for complex network characterization [49,64,65] are as follows: The k i of v i refers to the number of edges connected to the node in the network. In a directed network, the k i is composed of two components: the k out i , which is equal to the number of outgoing edges or successor nodes, and the k in i , which is equal to the number of incoming edges or predecessor nodes. The k i is the simplest and most effective concept for measuring the importance of a node in a network: where e ij is a connected edge from v i to v j . If this edge exists, e ij = 1; otherwise, e ij = 0. Similarly, e ji is a connected edge from v j to v i , and if this edge exists, then e ji = 1; otherwise, e ji = 0. P(k) is defined as the probability of selecting a node randomly in the network whose degree is k, i.e., the ratio of the number of nodes with k i = k to the number of all nodes in the network. P(k) describes the most basic topological characteristics of the network and is also an important scale used to identify network types: where n k is the number of nodes in the network with the degree k, and n is the order of the adjacency matrix (i.e., the number of nodes in the network). CC i represents the connection relationship among all nodes that are adjacent to a node in the network. CC i is defined as the proportion of the actual number of connected edges of these neighboring nodes to the maximum number of possible connected edges. CC i reflects the importance and subsequent growth of this node among neighboring nodes. CC is defined as the mean value of the CC i of all nodes in the network, which is an important index to measure network collectivization and represent the clustering ability of the network: where e i is the actual number of connected edges of neighboring nodes for v i , and the interconnection between two nodes in a directed network serves as two edges. m i is the number of neighboring nodes of v i , which is defined as the sum of successor and predecessor nodes in a directed network (repeated nodes are only counted once). BC i represents the probability of the shortest paths through this node in the network, and it is defined as the proportion of the number of shortest paths through this node to all shortest paths. BC i reflects the load and influence of a particular node in the network, and network robustness can be assessed by attacking nodes with high BC i : where n st is the number of shortest paths from v s to v t that pass through v i , and N st is the total number of shortest paths from v s to v t . The degrees of dispersion of BC i for different nodes are large, thus normalization processing is performed: l ij is defined as the sum of edge weights for the shortest path between v i and v j (a pair of nodes) in the network. l ij is used to measure the shortest distance between two nodes in a network, which, in this paper, represents the maximum connectivity between two nodes. The maximum value of l ij between any two nodes v i and v j is denoted as D of the whole network. L reflects the overall connectivity and connection efficiency of the network, which is of great significance. L is defined as the mean value of l ij between two random connected nodes in a directed network: where n l is the number of connected node pairs in the network.
To avoid divergence in the calculation of L caused by unconnected node pairs, the concept of global efficiency is proposed. GE is negatively correlated with L, which can quantify the efficiency of the network in transmitting information between nodes. GE is an important index to measure the connectivity of one network, and is defined as the mean value of the reciprocal of l ij between two random nodes in the network: If there is an interconnected path between v i and v j in the directed network, both nodes are strongly connected. If every two nodes in a network are strongly connected, this directed network graph is a strongly connected graph. The subgraph of a directed network graph with the maximum connection strength is called the strongly connected component, and the strongly connected component with the largest scale is called the maximum strongly connected component. MK is defined as the number of nodes in the maximum strongly connected component of the directed network graph, which reflects the robustness of the network under attack.
As the sum of edge weights is needed to calculate the shortest path, the concept of entropy is used to represent edge weights in the network, to meet the needs of studying the network characteristics and searching paths:

Shortest Paths of CN
There are two commonly used algorithms for shortest paths in CN: the Floyd algorithm and the Dijkstra algorithm. The Floyd algorithm [66] adopts the idea of dynamic program-ming, and it is suitable for solving the shortest path problem between any two nodes. The Dijkstra algorithm [67] adopts a greedy search strategy, and it is suitable for solving the shortest path problem with a single source. The Yen algorithm [68] is currently the most widely used algorithm to solve the problem of K shortest paths. Considering the lower time complexity of the Dijkstra algorithm, the shortest path is calculated by the Dijkstra algorithm. However, the conventional Dijkstra algorithm can only calculate one shortest path for the same initial node and target node. Therefore, the HADY is proposed to calculate all shortest paths of specified node pairs in the network.
The idea of HADY is firstly to obtain a shortest path of the specified node pair in the original network by the Dijkstra algorithm. By constantly removing the edges and their combinations of known shortest paths from the original network, the shortest path network of the specified node pair is obtained, and then, the shortest path network is transformed into an unweighted network. Finally, the Yen algorithm is used to solve K shortest paths of the shortest path network, and multiple shortest paths, sorted by the number of nodes, are obtained. The flow chart of HADY is shown in Figure 2.
the network characteristics and searching paths: There are two commonly used algorithms for shortest paths in CN: the Floyd algorithm and the Dijkstra algorithm. The Floyd algorithm [66] adopts the idea of dynamic programming, and it is suitable for solving the shortest path problem between any two nodes. The Dijkstra algorithm [67] adopts a greedy search strategy, and it is suitable for solving the shortest path problem with a single source. The Yen algorithm [68] is currently the most widely used algorithm to solve the problem of K shortest paths. Considering the lower time complexity of the Dijkstra algorithm, the shortest path is calculated by the Dijkstra algorithm. However, the conventional Dijkstra algorithm can only calculate one shortest path for the same initial node and target node. Therefore, the HADY is proposed to calculate all shortest paths of specified node pairs in the network.
The idea of HADY is firstly to obtain a shortest path of the specified node pair in the original network by the Dijkstra algorithm. By constantly removing the edges and their combinations of known shortest paths from the original network, the shortest path network of the specified node pair is obtained, and then, the shortest path network is transformed into an unweighted network. Finally, the Yen algorithm is used to solve K shortest paths of the shortest path network, and multiple shortest paths, sorted by the number of nodes, are obtained. The flow chart of HADY is shown in Figure 2.

SP=Ø; k=1;
The initial node is vs; The target node is vt  The detailed steps of this algorithm are as follows: 1. Set the edge set of shortest paths SP = Ø (the initial SP is the empty set). The initial node is v s , and the target node is v t . The initial k is 1. 2. The Dijkstra algorithm is used to obtain a shortest path from v s to v t in the original network G. All edges of the shortest path are put into the set SP. The number of edges is denoted as sn, and the length of the path is denoted as l. 3. List all combinations of k edges taken from the set SP and put them into the set SK.
There are kn combinations in total, and the initial i is 1. 4. Take the i-th combination from SK and delete the set of corresponding edges SK(i) in the original network G. The Dijkstra algorithm is used to find the shortest path from v s to v t in the modified network G 1 , and all edges of the shortest path are put into the set ST. The length of the path is denoted as d, and i = i + 1. 5. Determine whether a new edge of the shortest path is added. If d = l and SP ∪ ST = SP, a new edge has been added. Then, the logical value a = 1, SP = SP ∪ ST, and the number of edges is recalculated and denoted as sn. Otherwise, no new edge is added, the logical value a = 0, SP = SP, and the number of edges is denoted as sn.
6. Determine whether all combinations have been taken out. If i ≤ kn, not all combinations are taken out, and step 4 will be taken. Otherwise, if all combinations have been taken out, step 7 will be taken. 7. Determine whether this combination should be carried out again. If a = 1, this combination should be carried out again, then k = 1. Otherwise, there is no need to carry out this combination again, then k = k + 1. 8. Determine whether all combinations have been listed. If k ≤ sn, not all combinations have been listed, then step 3 will be taken. Otherwise, all combinations have been listed, then step 9 will be taken. 9. The set of all edges of shortest paths SP is used to form the shortest path network graph G 2 . All edge weights in G 2 are set to 1 (i.e., unweighted network), and K = 1. 10. K shortest paths from v s to v t in the unweighted shortest path network graph G 2 are obtained by the Yen algorithm and put into the set SP K . 11. Determine whether the combined algorithm has been completed. If SP K = Ø, step 12 will be taken; otherwise, K = K + 1 and step 10 will been taken. 12. The sets SP 1 to SP K are all the shortest paths from v s to v t in the original network G, and the length of the shortest path is l.

Case Study
A case study of major offshore oil and gas industry accidents caused by storm disasters is carried out, applying the proposed method. Firstly, causal factors are identified. Then, the hierarchy structure of the evolution of causes is studied. Finally, the evolution model of the accident causes is constructed to analyze the network characteristics and calculate the shortest path. The flow chart of the case study is shown in Figure 3. In this paper, a total of 78 major accidents in the offshore oil and gas industry caused The selection of original data is the first step of the research, and its accuracy and detail are the foundation of successful research, which also guarantees the reliability of the results. To guarantee the accuracy and reliability of the results, the data selected in this paper are mainly obtained from official accident investigation reports, accident announcements and accident statistics of competent authorities, news reports of mainstream media, and accident databases of authoritative institutions. Meanwhile, we only selected accidents with detailed processes, and we focused on those cases that cause serious structural damage or casualties. In addition, three groups of experts (from university, research institute, and offshore platform) were responsible for the analysis. Each team consisted of three members and performed an independent analysis of the same data. Through comparison and analysis, the final results were extracted to avoid inconsistency of concepts and categories caused by differences in theoretical sensitivity. Through these measures, the accuracy, reliability, and availability of the results are further guaranteed.

Coding of Accident Causal Factors
In this paper, a total of 78 major accidents in the offshore oil and gas industry caused by storms were selected, involving 215 different types of data. Details of the accidents are presented in Appendix A [69][70][71][72][73][74]. Among the selected 78 cases, 72 cases were randomly selected for coding, and the remaining six cases were used for theoretical saturation tests. After coding of the accident-causing factors based on the selected accident cases, the results of three groups of analysis were compared, discussed, and integrated. Subsequently, a total of 63 concepts, 25 categories, nine main categories, and three core categories were obtained, as shown in Table 1. Finally, an experienced researcher was invited to conduct independent coding of the remaining six cases (No. 5,No.12,No.18,No.31,No.42,and No.66) based on the same method. It was found that no new coding was generated, proving that the existing coding has covered all the accident cases. This indicated that the results satisfy the theoretical saturation test requirements. If the new coding appears, this indicates that the existing coding does not cover all accident cases. Then more accident cases need to be added and recoded until the verification test is passed. This step can further ensure the accuracy and reliability of the analysis results.  According to Table 1, inclement weather (y 1 ), loss of watertight integrity (y 13 ), lack of risk awareness (y 21 ), poor state of the platform (y 14 ), damage of support structures (y 10 ), damage of platform facilities (y 12 ), and damage of watertight structures (y 11 ) are the most frequently occurring categories. Attention needs to be focused in those categories.

Analyzing Interaction Relationships among Causal Factors
To analyze interaction relationships among different categories, the matrix S = [s ij ] n×n is defined, where s ij (i = 1, 2, . . . , n; j = 1, 2, . . . , n) represents the interaction value between the two categories. s ij is defined as the ratio of the number of intersection cases to the number of union cases the two categories belong to. After calculating and processing, the matrix R = [r ij ] n×n is obtained. {0, 1, 2, 3} in the matrix represents {no interaction, weak interaction, medium interaction, strong interaction}, respectively. The calculation rules are as follows:

Developing the Hierarchy of Causes of Offshore Storm Accidents
Based on the 25 causal factors (i.e., the categories) extracted by GT, DEMATLE-ISM is used to analyze influencing relationships and establish the hierarchical structure of the accident causation evolution.

Analyzing Influence Relationships among Causal Factors
The centrality and causality of accident-causing factors are plotted on a Cartesian coordinate system, as shown in Figure 4.  According to Figure 4, factors such as damage to watertight structure (y11), damage to platform facility (y12), loss of water-tightness (y13), and poor state of the platform (y14) have higher centrality. This indicates that these factors play an important role in the accident-causing process. Inclement weather (y1), lack of risk awareness (y21), and inadequate According to Figure 4, factors such as damage to watertight structure (y 11 ), damage to platform facility (y 12 ), loss of water-tightness (y 13 ), and poor state of the platform (y 14 ) have higher centrality. This indicates that these factors play an important role in the accidentcausing process. Inclement weather (y 1 ), lack of risk awareness (y 21 ), and inadequate safety training (y 22 ) have higher causality, indicating that these factors greatly impact other factors in the process of accident development and are the main causes of accidents.

Establishing the Hierarchical Structure of Causal Factors
In the transformation process from the matrix H to the matrix K, λ = 0.09 is chosen as the appropriate threshold. According to the obtained K, the system hierarchy is divided, and then, the accident evolution hierarchy is determined. The result is shown in Figure 5, in which only the influencing relationships among factors at adjacent levels are presented. In Figure 5, the circle shape represents the cause factor, while the hexagon shape represents the result factor. The size represents the category frequency obtained in the analysis of GT, which is divided into four levels (the larger the size, the higher the frequency of occurrence). The shading represents the centrality of the factor, which is divided into three levels (the darker the color, the greater the centrality). The thickness of the line between different levels represents the influence relationship between levels (the thicker the line, the stronger the influencing relationship). According to Figure 4, factors such as damage to watertight structure (y11), damag to platform facility (y12), loss of water-tightness (y13), and poor state of the platform (y14 have higher centrality. This indicates that these factors play an important role in the acci dent-causing process. Inclement weather (y1), lack of risk awareness (y21), and inadequat safety training (y22) have higher causality, indicating that these factors greatly impac other factors in the process of accident development and are the main causes of accidents

Establishing the Hierarchical Structure of Causal Factors
In the transformation process from the matrix H to the matrix K, λ = 0.09 is chosen a the appropriate threshold. According to the obtained K, the system hierarchy is divided and then, the accident evolution hierarchy is determined. The result is shown in Figure 5 in which only the influencing relationships among factors at adjacent levels are presented In Figure 5, the circle shape represents the cause factor, while the hexagon shape repre sents the result factor. The size represents the category frequency obtained in the analysi of GT, which is divided into four levels (the larger the size, the higher the frequency o occurrence). The shading represents the centrality of the factor, which is divided into thre levels (the darker the color, the greater the centrality). The thickness of the line between different levels represents the influence relationship between levels (the thicker the line the stronger the influencing relationship). According to Figure 5, the 25 causal factors are divided into six layers. L1 is the sur face layer, which is the direct cause of the accident. L2 and L3 form the shallow layer, and represent early signs of the accident. L4 and L5 form the deep layer, and represent th According to Figure 5, the 25 causal factors are divided into six layers. L1 is the surface layer, which is the direct cause of the accident. L2 and L3 form the shallow layer, and represent early signs of the accident. L4 and L5 form the deep layer, and represent the concentrated emergence of accident causes. L6 is the ground layer, which is the root of the accident. In terms of factors, inclement weather (y 1 ) and lack of risk awareness (y 21 ) are the root causes of accidents. Collision of fixed object at sea (y 4 ), collision of floating object at sea (y 5 ), and collision of object on the platform (y 6 ) have complicated influence relationships with other factors that belong to upper and lower levels. Damage to support structure (y 10 ) and loss of water-tightness (y 13 ) have high centrality and occurrence frequency, and are key factors in accidents. Damage to platform facility (y 12 ) is the only path for the propagation of the causative factors at upper and lower layers. When the watertight structure is damaged (y 11 ) and the platform is in poor condition (y 14 ), accidents are more likely to happen than usual. At this time, all necessary control measures should be taken, and emergency work should be carried out to mitigate the consequences of the accident.

Analyzing the Evolution Process of Offshore Storm Accidents
Based on the modeling principle of CN, the accident evolution model of the offshore oil and gas industry under storm disasters is established, combining 63 concepts extracted from GT by analyzing the evolution hierarchy of accidents. Using this model, the accident evolution characteristics are analyzed, and the shortest evolution paths are calculated.

Modeling the Evolution Process of Accident Causes
In line with the construction method of the matrix R, the matrix RX is established. In combination with the matrix F, 63 concepts are put into the set of X (x i &x j ∈ X) to establish the matrix A = [a ij ] n×n , where n represents the order of A, and a ij represents the connectivity from x i to x j . In the matrix A, {0, 1, 2, 3} represents {no connectivity, weak connectivity, medium connectivity, strong connectivity}, respectively. According to the matrix A, a weighted directed CN graph is plotted, in which v i is represented by x i , and a ij is taken as w ij . In the established graph, the thickness and length of the line represent the connectivity between two nodes. Corresponding to thicker and shorter lines, a larger weight represents a greater connectivity, as shown in Figure 6. emergency work should be carried out to mitigate the consequences of the accident.

Analyzing the Evolution Process of Offshore Storm Accidents
Based on the modeling principle of CN, the accident evolution model of the offshore oil and gas industry under storm disasters is established, combining 63 concepts extracted from GT by analyzing the evolution hierarchy of accidents. Using this model, the accident evolution characteristics are analyzed, and the shortest evolution paths are calculated.

Modeling the Evolution Process of Accident Causes
In line with the construction method of the matrix R, the matrix RX is established. In combination with the matrix F, 63 concepts are put into the set of X (xi&xj ∈ X) to establish the matrix A = [aij]n×n, where n represents the order of A, and aij represents the connectivity from xi to xj. In the matrix A, {0, 1, 2, 3} represents {no connectivity, weak connectivity, medium connectivity, strong connectivity}, respectively. According to the matrix A, a weighted directed CN graph is plotted, in which vi is represented by xi, and aij is taken as wij. In the established graph, the thickness and length of the line represent the connectivity between two nodes. Corresponding to thicker and shorter lines, a larger weight represents a greater connectivity, as shown in Figure 6. In Figure 6, the whole network is obviously divided into two subnetworks, which are the region A centered on x1 (strong ocean storm) and huge x2 (ocean wave), and the region B centered on x51 (inadequate risk perception) and x52 (inadequate risk response). The internal connection of these two subnetworks is complex and the average connectivity is high, which means that internal factors are closely related.  In Figure 6, the whole network is obviously divided into two subnetworks, which are the region A centered on x 1 (strong ocean storm) and huge x 2 (ocean wave), and the region B centered on x 51 (inadequate risk perception) and x 52 (inadequate risk response). The internal connection of these two subnetworks is complex and the average connectivity is high, which means that internal factors are closely related. The region A and the region B are connected through intermediate nodes including lack of x 17 (watertight subdivision isolation), x 47 (violation of towing operation), et al. The nodes at the edge of the subnetwork have higher connectivity with intermediate nodes, but the internal nodes of the two subnetworks need to be connected through edge nodes as medium. The average connected path between these nodes is long and the connection relationship between them is sparse. The whole network shows clear community structure characteristics.

Analyzing the Evolution Characteristic of Accident Causes
The k i , k in i , and k out i in the network are calculated, and nodes with high values are: x 1 (strong ocean storm), x 2 (huge ocean wave), x 8 (collision of buildings), x 9 (collision of vessels), x 12 (collision of cargos), x 13 (collision of lifeboats), x 37 (list of the platform) and x 51 (inadequate risk perception). At the same time, P(k) in the network is calculated, which presents the characteristics of power law distribution (P(k)~k −β ) where β = 2.82, which is consistent with the distribution feature of β = 2~3 in a scale-free network. The CC i and BC i of all nodes in the network are calculated, as shown in Table 2. According to Table 2 and based on Figure 6, the network characteristic graph of the accident evolution is drawn, as shown in Figure 7. The size of the node represents the value of CC i , and the color shade of the node represents the value of BC i . A larger size represents greater CC i , and a darker shade represents greater BC i . consistent with the distribution feature of β = 2~3 in a scale-free network. The CCi and BC of all nodes in the network are calculated, as shown in Table 2.  Table 2 and based on Figure 6, the network characteristic graph of the accident evolution is drawn, as shown in Figure 7. The size of the node represents the value of CCi, and the color shade of the node represents the value of   In Table 2 and Figure 7, the maximum CC i is 0.5, and the corresponding nodes are x 15 (improper screws), x 19 (design defect of single hull tankers), x 44 (poor performance of protective clothing), and x 48 (wrong wellhead connection). However, the number of neighboring nodes of these nodes is too small to be representative. Among the other nodes, the maximum CC i is 0.45, corresponding to x 55 (insufficient auxiliary vessels), x 56 (insufficient lifeboats), and x 57 (insufficient emergency protective equipment), and these nodes have many neighboring nodes, indicating their relatively high collectivization. Therefore, adequate emergency equipment can effectively suspend the rapid escalation of accidents. The CC of the network is 0.223, which is far less than 1, but also much more than 1/n = 0.016. This is in line with the characteristics of the real-world network, which proves the feasibility of the established network.
In Table 2 and Figure 7, nodes with high BC i are x 41 (power system failure), x 38 (thruster failure), x 37 (list of the platform), x 31 (flooding of deck), and x 33 (flooding of subdivision), indicating that there are multiple shortest paths through these nodes. In other words, these nodes are where the accident evolution most likely passes. Therefore, the proper functioning of the power system, thrusters, and other equipment should be ensured, as this can effectively reduce the network connectivity and the possibility of accidents. The same applies to the stability and water tightness of the platform.
The L of the network is 12.69, and the D of the network is 32.24. L is far less than D, indicating that the accident evolution has small-world network characteristics. The GE is 4.01 × 10 −2 , which indicates that the network has strong connectivity ability and efficiency. The MK is 21, indicating that interconnected paths exist among these 21 nodes in the network, and the overall robustness of the network is strong.

Calculating Shortest Evolution Paths of Accident Causations
According to the evolution hierarchical structure and statistical frequency shown in Figure 5 and Table 1, x 1 (strong ocean storm), x 2 (huge ocean wave), x 3 (rapid ocean current), x 51 (inadequate risk perception), and x 52 (inadequate risk response) in the ground layer (L6) are selected as initial nodes, and x 31 (flooding of decks), x 33 (flooding of cabins), and x 37 (list of the platform) in the surface layer (L1) are selected as target nodes. These nodes result in a total of 15 node pairs (i.e., 15 shortest path groups). The proposed HADY is used to obtain the shortest path network graphs, and the thickness of the line represents the connectivity between two nodes, where stronger connectivity corresponds to thicker lines, as shown in Figure 8. The network graphs of shortest paths for different target nodes are presented in Appendix B.
bility of the established network.
In Table 2 and Figure 7, nodes with high BCi are x41 (power system failure), x3 (thruster failure), x37 (list of the platform), x31 (flooding of deck), and x33 (flooding of sub division), indicating that there are multiple shortest paths through these nodes. In other words, these nodes are where the accident evolution most likely passes. Therefore, the proper functioning of the power system, thrusters, and other equipment should be en sured, as this can effectively reduce the network connectivity and the possibility of acci dents. The same applies to the stability and water tightness of the platform.
The L of the network is 12.69, and the D of the network is 32.24. L is far less than D indicating that the accident evolution has small-world network characteristics. The GE is 4.01 × 10 −2 , which indicates that the network has strong connectivity ability and efficiency The MK is 21, indicating that interconnected paths exist among these 21 nodes in the net work, and the overall robustness of the network is strong.

Calculating Shortest Evolution Paths of Accident Causations
According to the evolution hierarchical structure and statistical frequency shown in Figure 5 and Table 1, x1 (strong ocean storm), x2 (huge ocean wave), x3 (rapid ocean cur rent), x51 (inadequate risk perception), and x52 (inadequate risk response) in the ground layer (L6) are selected as initial nodes, and x31 (flooding of decks), x33 (flooding of cabins) and x37 (list of the platform) in the surface layer (L1) are selected as target nodes. These nodes result in a total of 15 node pairs (i.e., 15 shortest path groups). The proposed HADY is used to obtain the shortest path network graphs, and the thickness of the line represents the connectivity between two nodes, where stronger connectivity corresponds to thicker lines, as shown in Figure 8. The network graphs of shortest paths for different target nodes are presented in Appendix B.  According to the HADY, the shortest paths and lengths are calculated, and a total of 43 shortest paths are obtained. The specific results are shown in Table 3.
In Table 3, a total of five pairs of nodes have the least number of paths (only one path), and the node pair with the largest number of paths is path_3, which has eight shortest paths. For the same node pair, the more paths there are, the more difficult it is to control the evolution process. Because even if one of the paths (i.e., the causative chain) is blocked, there may be other evolution paths. Therefore, attention should be paid to node pairs with multiple shortest paths, and the node or edge with the highest blocking efficiency should be selected as control. For example, if the edge x 37 → x 31 or the node x 37 is removed from the eight shortest paths of path_3, the number of shortest paths decreases from eight to two, and the blocking efficiency increases to 75%.  42 43 path_15A path_15B In addition, the node pair of the minimum shortest path length is (x 2 , x 31 ) with a length of 2.30, while the node pair of the maximum shortest path length is (x 3 , x 31 ) with a length of 13.82. The shorter path length represents stronger connectivity and greater possibility of accident occurrence. All shortest paths have at least one edge and five edges at most. For the same path length, more edges represent more nodes, i.e., the development process of accidents involves more causative factors. This also indicates that there are more optional control measures.
In general, the longer the shortest path length of a node pair, the more paths and the more edges there are. However, high risk features (e.g., shorter path length, fewer edges, and more paths) generally do not appear on the same path. Therefore, when controlling the evolution path of accidents, the high risk features of paths should be identified first. Subsequently, appropriate prevention and control measures should be taken.
To quantify the importance of different shortest path groups, according to the 15 shortest path groups listed in Table 3, intentional attacks of two modes are carried out on the original network: 1. Attack only the edges of the shortest path; 2. Attack both the edges and intermediate nodes of the shortest path.
The l, GE and MK under the two attack modes are analyzed, as shown in Figure 9.
To quantify the importance of different shortest path groups, according to the 15 shortest path groups listed in Table 3, intentional attacks of two modes are carried out on the original network: 1. Attack only the edges of the shortest path; 2. Attack both the edges and intermediate nodes of the shortest path.
The l, GE and MK under the two attack modes are analyzed, as shown in Figure 9. In Figure 9a, the l of different shortest path groups under the attack mode 1 is consistent with that under the attack mode 2. Compared with the original network, all l increase under both attack modes. The average increased proportion of all l is 35.19%, indicating that the shortest path length increases and the connectivity reduces under both attack modes, which can reduce the possibility of accidents. In Figure 9b In Figure 9a, the l of different shortest path groups under the attack mode 1 is consistent with that under the attack mode 2. Compared with the original network, all l increase under both attack modes. The average increased proportion of all l is 35.19%, indicating that the shortest path length increases and the connectivity reduces under both attack modes, which can reduce the possibility of accidents. In Figure 9b, the average decreased proportion of GE is 3.43% under the attack mode 1, and the average decreased proportion of GE is 8.94% under the attack mode 2. A comparison of both attack modes shows that the change trend of GE is almost identical. In Figure 9c, the average MK of the network is 20.73 under the attack mode 1, and the average MK of the network is 18.1 under the attack mode 2. Compared with the attack mode 1, the MK of the network decreases most when attacking path_3 under the attack mode 2, and the decrease proportion is 42.11%. Compared with the original network, the attack effects of the two attack modes are not ideal. This indicates that the original network has strong robustness. In such a case, attacking a single path or node can hardly reduce the robustness of the network.
The research above treats multiple shortest paths of the same node pair as a whole. To study the influences of different shortest paths of the same node pair on the network, it is necessary to integrate various network characteristics under attack. After normalization and summation, the importance of different shortest paths PS for the same node pair is obtained. The calculation formula is as follows: where GE and MK are the values after attacking the shortest path. MS is the importance of intermediate nodes for the shortest path. GE max , MK max , and MS max are the maximum values of all shortest paths for the same node pair. GE min , MK min , and MS min are the minimum values of all shortest paths for the same node pair. According to Equation (17), the importance of different shortest paths for the same node pair are calculated and sorted in the same node pair, as shown in Table 4.

Results of the Case Study
According to the above analysis results, the shortest path groups that need to be focused on are as follows: 1. The path_3 (x 3 , x 31 ), which is the shortest path group hardest to prevent, i.e., the one with the most paths. There are eight shortest paths, and each path contains the initial node x 3 (rapid ocean current) and the target node x 31 (flooding of decks). The intermediate nodes are x 9 (collision of vessels), x 10 (collision of debris), x 11 (collision of cantilever deck), x 23 (fracture of legs), x 24 (collision of vessels) and x 37 (list of the platform). 2. The path_2 (x 2 , x 31 ), which is the shortest path group most likely to occur, i.e., the one with the minimum shortest path length. The path length is 2.30 and the initial node x 2 (huge ocean wave) of this path is directly connected to the target node x 31 (flooding of decks) without passing through any intermediate nodes. 3. The path_5 (x 52 , x 31 ) and path_14 (x 51 , x 37 ), which are the shortest path groups hardest to block, i.e., the one with the most edges. There are four shortest paths with five edges: the path_5B, the path_5C, the path_14F and the path_14G.
According to the attack effect analysis on different shortest path groups, the attack mode 2 achieves the better effects. In other words, to prevent the occurrence of accidents, it is necessary not only to control the occurrence of risk events (i.e., nodes in the network), but also to block the propagation path between risk events. In addition, due to the limitations imposed by cost, it is impossible to guard against all risk events in the production and operation activities. For the same amount of time and manpower spent, safety investments should be skewed towards risky events that are likely to yield higher safety benefits.
Based on the analysis results of l, GE and MK after attacking different node pairs, the risky events in the following shortest path groups deserve more safety investments: (a) path_3, which consists of eight paths; (b) path_5, which consists of three paths; (c) path_14, which consists of seven paths; (d) path_6, which cconsists of four paths; (e) path_15, which consists of two paths.
In line with the results in Table 4, when some shortest paths are attacked, the connectivity and robustness of the network are greatly reduced. Meanwhile, these paths can also interfere with other shortest paths, indicating that they are critical paths of the entire network. In conclusion, the focus should be on preventing and controlling the following risk events and propagation paths: 1. The path_3E, x 3 (rapid ocean current) → x 9 (collision of vessels) → x 23  In the above five shortest paths, three paths all contain x 52 → x 53 → x 17 → x 33 → x 37 . It is proved that this shortest path plays an important role in the whole accident evolution network. Therefore, controlling the propagation of this path has great significance for preventing the major storm-disaster-induced accidents in the offshore oil and gas industry.
According to the above analysis results of accident evolution path, different measures are proposed for different stages of accidents. Specific measures are as follows: 1. Preventive measures before accidents: accurate storm forecast and early warning; redundant design of structure and equipment; regular detection of structural defects; fixing movable items on the deck; ensuring the strength of towing ropes; ensuring water-tightness; enhancing risk awareness and strengthening safety training; complete operation regulations and emergency plans; adequate emergency equipment. 2. Control measures in accidents: closing sea valves on the platform; timely drainage of cabin and deck; adjusting ballast conditions of the platform; avoiding excessive oscillation and list of the platform; avoiding hot work; disconnecting from the wellhead; ensuring the normal operation of the power system; ensuring the normal operation of the dynamic positioning system; timely strengthening the support structure; stay out of the storm zone by towing or self-propulsion. 3. Emergency measures after getting out of control: activating the alarm system; emergency muster at the designated location; evacuating platform personnel in sequence; calling for outside help; shutting down unessential electrical equipment; distributing personal survival equipment; lowering lifeboats, life rafts and emergency escape ladders; rescue assist of guard ships.

Discussion
Compared with the existing literature, on the one hand, the introduction of GT compensates for the shortcomings of the traditional identification of accident causes, which relies too much on the knowledge and experience levels of personnel. GT is characterized by strong objectivity, which can ensure the accuracy and consistency of classification results of accident cause factors. On the other hand, most of the current methods usually divide the causal factors into large categories first and then determine small categories, which leads to the unstable level and content of classification. According to programmatic GT in this paper, a total of 63 concepts, 25 categories, nine main categories, and three core categories were obtained in turn. The result shows that this inductive method can solve the problem well.
Meanwhile, the method of DEMATEL-ISM is improved in this paper. In the establishment of the direct influence matrix, the interaction matrix obtained from the objective analysis in the previous step is introduced based on expert rating. As a result, the strong subjectivity associated with the traditional establishment of the matrix can be avoided. The combination of these two methods retains the advantages of objective induction and subjective reasoning, respectively.
At the same time, the analysis results of the first two steps provide input for the construction of the CN model, and the three methods form a step-by-step systematic evaluation method. In addition, a hybrid algorithm named HADY is proposed. In previous studies, when calculating the shortest path of CN model, each node pair usually only gets one path. However, the same node pair may have multiple shortest paths. Using HADY, multiple shortest paths under the same node pair can be calculated effectively. We obtained seven shortest paths on the node pair path 14. This is proved that the algorithm is effective and more in line with the actual situation.
Through the actual case study of the proposed method, the importance of different nodes and paths is calculated, and five key nodes and five critical paths are obtained. They are the node x 1 (strong ocean storm), the node x 9 (collision of vessels), the node x 31 (flooding of decks), the node x 37 (list of the platform), the node x 51 (inadequate risk perception), the path_3E (x 3 → x 9 → x 23 → x 37 → x 31 ), the path_5C (x 52 → x 53 → x 17 → x 33 → x 37 → x 31 ), the path_14G (x 51 → x 52 → x 53 → x 17 → x 33 → x 37 ), the path_6C (x 1 → x 9 → x 26 → x 33 ), and the path_15B (x 52 → x 53 → x 17 → x 33 → x 37 ). Compared with the research results in the existing literature, the proposed method can not only obtain the importance of nodes, but also obtain the importance of different paths through the method of path attack. The results provide more references for accident prevention.
In general, the proposed method contributes to the precise identification of accident causes and the systematic analysis of accident development. Based on the results of the analysis, reasonable and effective control measures can be implemented to interfere with the development process of accidents. This will slow down or even control the evolution process of accidents under storm disasters to the greatest extent and ensure the safety of offshore oil and gas operations. Further, the application of research results will help reduce casualties caused by major accidents. Meanwhile, it is helpful to achieve sustainable development by reducing the impact of oil and gas accidents on the marine environment.
However, the method proposed in this paper also has its limitations. For example, the method requires detailed data and takes a long time to prepare and code. Moreover, the current CN model of accident evolution cannot realize the dynamic adjustment of nodes and connecting edges. Further research can integrate artificial intelligence in the extraction process of accident causal factors to realize intelligent coding of causal factors for many accident cases. How to embody the dynamic changes of CN models is also one of the important tasks of our future research.

Conclusions
This paper presents a method for studying the evolution process of accident causes by combining GT, DEMATEL-ISM, and CN. The integrated method covers different stages, including the identification of accident causes, their qualitative analysis, and quantitative assessment, which can be applied for systematically studying the development process of accidents. The application of the proposed method is illustrated by a case study. The results show that inclement weather and lack of risk awareness are the root causes of accidents. Moreover, the damage of support structures and the loss of water-tightness are key nodes of the accident development. Furthermore, the damage of platform facilities is the inevitable path of accident evolution. The constructed accident evolution network model conforms to the characteristics of small-world network and scale-free network. The overall connectivity and robustness of the network are high, and it is difficult to reduce the robustness of the network by attacking a single path or node. Results show that there are five key nodes and five critical paths in the process of accident evolution. Through targeted prevention and control of these nodes and paths, the average shortest path length of the accident evolution network is increased by 35.19%, and the maximum global efficiency decreases by 20.12%. Thus, more attention should be focused on the paths with higher importance as identified by the evaluation results of the attack effect under the limitations imposed by cost, time, and manpower.

Annotations
X-the set of concepts; Y-the set of categories; F-the direct influence matrix; C-the normal matrix; T-the comprehensive influence matrix; I-the identity matrix; H-the overall influence matrix; K-the reachable matrix; λ-the threshold value; b i -the influencing degree; d i -the influenced degree; c i -the centrality; a i -the causality; Q i -the antecedent set; P i -the reachable set; m i -the factor degree; S-the direct interaction matrix; R-the interaction matrix of categories; RX-the interaction matrix of concepts; A-the adjacency matrix; G-the CN graph model; V-the set of nodes; E-the set of edges; W-the set of weights; k i -the node degree; k out i -the out-degree; k in i -the in-degree; P(k)-the degree distribution; CC i -the clustering coefficient; CC-the network clustering coefficient; BC i -the betweenness centrality; BC i -the normalized betweenness centrality; l ij -the shortest path length; D-the network diameter; L-the average shortest path length; GE-the global efficiency; MK-the maximum strongly connected component scale; ew ij -the entropy weight.