TY - JOUR AU1 - Li, Qiaochu AU2 - Zhang, Peng AB - 1. Introduction China’s goals of carbon peaking and carbon neutrality have clarified the strategic direction and requirements for comprehensive green transformation [1]. As a clean, efficient and high-quality fossil energy, natural gas has become an important transitional energy source for the third energy transition, which drives the continuous construction of gas pipeline network. The karst landforms in China are widely distributed, concentrated in the southwest carbonate geological areas such as Sichuan, Guizhou, Yunnan, Hunan, Guangxi, and Hubei. The southwest China is rich in natural gas resources, and its underground gas pipeline network is widely distributed. In this regard, the gas pipelines inevitably pass through the karst areas, and the geological development poses significant hidden dangers to the safe and stable operation of buried pipelines [2]. Meanwhile, the natural environment in southwest China is complex and fragile. The terrain is mainly mountainous and plateau, with broken surface rocks and rainy climate. In addition, the gas pipeline disaster in karst areas is a special case of disaster. Specifically, it is dominated by natural factors, but coupled with human factors to form a more complex system structure. In this regard, karst collapse is easy to form under the coupling influence of natural events such as rainstorm and human activities such as third-party construction. It will lead to the failure and rupture of buried gas pipelines, and further induces secondary disasters such as fires and explosions. Those disasters have a negative impact on energy supply, social stability and economic development, and seriously affect the high-quality transformation of China’s energy industry [3]. Therefore, it is necessary to scientifically and accurately evaluate the disaster-causing probability of gas pipelines in karst areas, and provide theoretical guidance for disaster emergency response and rescue work. The evaluation of disaster-causing possibility includes failure analysis and probability calculation [4]. In recent years, pipeline disasters caused by soil collapse have occurred frequently. Scholars have begun to pay attention to the failure analysis of buried pipelines in soil collapse areas, and computer numerical simulation methods with advantages such as high accuracy and strong realism have been widely used [5]. The numerical simulation method sets model parameters and boundary conditions based on the actual engineering characteristics, and further establishes a three-dimensional finite element calculation model by combining intelligence analysis software such as ABAQUS and ADINA. In this regard, the mechanical response of pipelines under different collapse scenarios and different operating conditions can be simulated, so as to accurately investigate the failure characteristics and disaster-causing laws of buried pipelines in soil collapse areas [6]. Trifonov [7] used the finite element analysis method to evaluate the stress-strain response of the cross-fault pipelines in the state of uniform internal pressure based on two different types of fault representations. Ukan et al. [8] studied the mechanical response of cross-fault pipelines with different diameter-thickness ratios and steel grades. Naeen et al. [9] investigated the deformational behavior of buried pipelines under soil collapse by applying 3D continuum finite element modeling. With the development of disaster probability-causing evolution for pipelines in soil collapse areas, in addition to technical analysis of pipeline failure, attention has been paid to the calculation of disaster-causing probability combined with multiple factors such as society, economy, and environment [10]. The commonly used models include fault tree [11], Bowtie node [12] and Bayesian network [13]. Meanwhile, on the basis of system structure analysis, the analytic methods such as cluster analysis, analytic hierarchy process, fuzzy analysis and system analysis are further used to carry out probability quantitative calculation [14]. Phan et al. [15] established a strain failure probability calculation model for buried pipelines crossing strike-slip faults based on the artificial neural network method. However, the predicted strain may be significantly different from the actual strain, thus significantly affecting the accuracy of possibility evolution. Aslkhalili et al. [16] established the calculation method of failure probability of buried pipelines under the influence of large lateral soil displacements based on the first-order reliability method, and used Python programming software to develop calculation code. Alvarado-Franco et al. [17] established the failure probability analysis model of buried pipelines under the influence of landslide based on the Monte Carlo simulation method, and conducted a case study using a pipeline in central Colombia. Zheng et al. [18] established the failure probability calculation model of buried pipelines under the influence of permanent ground displacement based on the finite difference method. However, the strain capacity was determined using existing equations in the literature. Compared with the numerical simulation results, its accuracy is insufficient. The above studies provides an effective reference for the disaster-causing probability evaluation of gas pipelines in soil collapse areas, but they still have limitations in the following aspects: firstly, the existing studies mainly focus on the mechanical failure process of pipelines under the influence of soil displacement, and conducts quantitative analysis based on engineering parameters such as collapse size, pipeline operation characteristics, and site soil properties [9, 19, 20]. The comprehensive analysis combined with disaster system theory is relatively lacking, and multiple factors such as society, economy and environment are not comprehensively considered. Therefore, it is often difficult to achieve accurate measurement of disaster-causing probability. Secondly, most of the existing studies only focus on the hazard analysis of disaster events on buried pipelines, without considering the pipeline vulnerability [21–23]. In this regard, the important influence of pipeline anti-external load ability on disaster-causing process is ignored. In addition, the collapse scenario mainly focuses on earthquake collapse, landslide, mining collapse, etc. It is worth noting that karst collapse has the characteristics of hiddenness and suddenness, and once it happens, it will have a devastating impact on the pipeline system. Therefore, the targeted research on it needs to be further strengthened. In this regard, this study proposes a method to analyze the disaster-causing probability of gas pipelines in karst areas based on disaster system theory and vulnerability theory. On the one hand, the hazard evaluation index system of gas pipeline disaster events in karst areas is established from multiple dimensions including the activity of disaster-prone environment, the risk of disaster factors and the vulnerability of disaster-bearing bodies. Combined with the advantages of information transmission and updating of the Bayesian network model, the hazard probability of disaster events is gradually calculated from a multi-level perspective. On the other hand, considering the pipeline vulnerability level and its relationship with the disaster type and intensity, pipeline structure and function, disaster and pipeline space-time configuration, etc., a mathematical model system is established based on the finite element simulation results and polyethylene (PE) pipeline characteristics. Combined with the design checkpoint method, the vulnerability probability of gas pipelines is calculated from the perspective of the interaction between disaster hazard intensity and pipeline resistance capacity. Finally, a new disaster-causing probability evaluation approach for gas pipelines in karst areas is formed, which comprehensively considers the event hazard and the pipeline vulnerability. It is helpful to expand the research ideas in the field of risk assessment of gas pipeline disasters, and provide theoretical basis for scientific, accurate and comprehensive investigation of disaster-causing probability of gas pipelines in karst sinkhole prone areas. 2. Methods 2.1. Hazard analysis of disaster events The hazard probability of disaster events is defined as the occurrence probability of all types of disaster events that threaten the safety of human beings and material wealth in the area. That is, the probability distribution of multiple risk inducing factors turning into serious disaster events [24]. 2.1.1. Establishment of hazard evaluation index system. Based on the disaster system theory, the hazard evaluation index system is established from three aspects: the activity of disaster-prone environment, the risk of disaster factors and the vulnerability of disaster-bearing bodies. (1) The disaster-prone environment. In a broad sense, the disaster-prone environment refers to a comprehensive social-natural system composed of multiple factors of atmosphere, lithosphere, biosphere, hydrosphere, and social material and cultural sphere. However, the formation of the disaster-prone environment is not simply the direct superposition of the related factors, but depends on the material circulation, energy flow and information transmission among various factors. The disaster-prone environment of gas pipeline disasters in karst areas is the breeding place of disaster factors, and it is also the indirect medium of the interaction between disaster factors and disaster-bearing bodies. The activity of disaster-prone environment (a) is closely related to the personnel unsafe state (k1), pipeline unsafe state (k2), environmental unsafe state (k3) and management loopholes (k4), as shown in Fig 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. The structure of disaster-prone environment. https://doi.org/10.1371/journal.pone.0316820.g001 Specifically, the personnel unsafe state (k1) mainly includes the low level of education (t1), the lack of professional skills (t2), the weak awareness of pipeline protection (t3), weak safety awareness (t4) and insufficient emergency response capacity (t5). The pipeline unsafe state (k2) mainly includes the pipeline aging (t6), pipeline corrosion (t7), the failure of safety facilities (t8) and construction technology defect (t9). The environmental unsafe state (k3) mainly includes the intensity of human activity (t10), construction and vibration (t11), groundwater extraction (t12), economic development level (t13), legal environment (t14), karst geological development (t15), groundwater activity (t16), overburden characteristics (t17), structural condition (t18), topographic features (t19) and meteorological condition (t20). The management loopholes (k4) mainly includes the lack of safety supervision (t21), the lack of safety publicity (t22), inadequate emergency support (t23) and unreasonable rules and regulations (t24). (2) The disaster factors. The disaster factors are the driving force and direct cause of hazards, which are manifested as various sudden extreme events in natural or social environments that threaten human life, cause material property losses, and cause ecological damage. They include both natural anomalies such as mudslides, tsunamis, tornadoes, and earthquakes, as well as anthropogenic anomalies such as desertification, hazardous spills, traffic accidents, and explosions. The disaster factors of gas pipelines in karst areas are the direct causes of casualties, economic losses and environmental damage, mainly including fire thermal radiation and explosion shock wave overpressure. The risk of disaster factors is mainly influenced by the probability, intensity, scope and duration of fire and explosion hazards. Among the various disaster factors caused by gas pipeline failure and leakage, the occurrence of a single disaster can lead to the subsequent occurrence of other secondary disasters. Therefore, the disaster factors of gas pipelines in karst areas show a chain formation process. The risk of disaster factors (b) is closely related to the fire risk (k5) and explosion risk (k6), as shown in Fig 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. The structure of disaster factors. https://doi.org/10.1371/journal.pone.0316820.g002 Specifically, the fire risk (k6) mainly includes the source of fire (t25), thermal radiation flux (t26), smoke and dust (t27), noxious fumes (t28), the types of gas leakage (t29), the scope of fire impact (t30), radiation time (t31), distance (t32) and space-limited situation (t36). The explosion risk (k5) mainly includes the shock wave overpressure (t33), noise (t34), the source of fire (t25), the types of gas leakage (t29), distance (t32), the scope of explosion impact (t35) and space-limited situation (t36). (3) The disaster-bearing body. The disaster-bearing bodies refers to the direct natural and social subjects affected by disasters, including various social resources such as human life, material property, social production and operation systems, as well as various natural resources such as animals, plants, soil, water sources, and atmosphere. The disaster-bearing bodies of gas pipeline disasters in karst areas are the direct object that bears the hazard of disaster factors, which mainly includes social, economic and natural disaster-bearing bodies along the gas pipelines. Different disaster-bearing bodies are closely related and interact with each other. The vulnerability of disaster-bearing bodies is mainly used to measure the difficulty of their being damaged by disasters. On the one hand, it depends on the physical characteristics of the disaster-bearing bodies, such as the strength and ductility level of gas pipelines. On the other hand, it depends on the cultivated disaster-bearing ability, that is, the multiple abilities of resistance, rescue and recovery in the face of disasters. The vulnerability of disaster-bearing bodies (c) is closely related to the social disaster-bearing body (k7), economic disaster-bearing body (k8), natural disaster-bearing body (k9) and disaster bearing capacity (k10), as shown in Fig 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. The structure of disaster-bearing bodies. https://doi.org/10.1371/journal.pone.0316820.g003 Specifically, the social disaster-bearing body (k7) mainly includes the human casualties (t37), the destruction of buildings and structures (t38) and the destruction of lifeline engineering (t39). The economic disaster-bearing body (k8) mainly includes the material property loss (t40), economic compensation (t41), resource waste (t42) and production and business suspension (t43). The natural disaster-bearing body (k9) mainly includes the atmospheric pollution (t44), soil hardening (t45), the damage to vegetation (t46) and the harm to animals (t47). The disaster bearing capacity (k10) mainly includes the early warning capability (t48), disaster resistance capability (t49), disaster relief capacity (t50) and recovery capacity (t51). 2.1.2. Establishment of Bayesian network model. The Bayesian network is a directed acyclic graph based on graph theory and probability theory, which is mainly used to describe the dependency relationships and probability distributions of a group of random variables. It consists of nodes, directed edges and probabilities, and a simple Bayesian network structure is shown in Fig 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. A simple Bayesian network structure. https://doi.org/10.1371/journal.pone.0316820.g004 In Fig 4, nodes A, B1, B2, C1 and C2 represent a group of random variables, and the arrows indicate their dependency relationships. The C1 points to B1, indicating that the C1 is the parent node of the B1, and the B1 is the child node of the C1. Meanwhile, the C1 and C2 have no parent node, so they are root nodes. The root nodes have corresponding prior probability distribution, and non-root nodes have corresponding conditional probability distribution. The Bayesian network can update the failure probability of the system by adjusting the probability distribution of nodes. Taking Fig 4 as an example, the adjustment process of probability distribution of nodes in the Bayesian network is investigated: assuming that at a certain moment t, the probability that the occurrence of C1 causes the occurrence of B1 is 0.2, the probability that the occurrence of C2 causes the occurrence of B1 is 0.3. The probability of B1 occurring when C1 and C2 both occur is 0.5, and the probability of B1 occurring when neither C1 and C2 occurs is 0.1. In this regard, the conditional probability distribution corresponding to the node B1 is shown in Table 1, and 0 and 1 respectively indicate the occurrence and non-occurrence of the root nodes. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. The conditional probability distribution corresponding to the node B1. https://doi.org/10.1371/journal.pone.0316820.t001 In the Bayesian network model, the probability calculation of event occurrence can be equivalent to the propagation and updating of beliefs [25]. In this regard, the calculation process is simple and the results are accurate. Therefore, the Bayesian network model can be used to evaluate the hazard probability of disaster events of gas pipelines in karst areas. Meanwhile, the fault tree model and Bayesian network model have similarities in principle, and the deductive reasoning function of the fault tree model can quickly determine the dependency relationships between multiple events. In this regard, the Bayesian network model can be established based on the corresponding mapping relationship. The transformation of the “or” gate and “and” gate in the fault tree model to the Bayesian network model is shown in Fig 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Schematic diagram of transformation between the fault tree model and the Bayesian network model. The “or” gate (A), and“and” gate (B). https://doi.org/10.1371/journal.pone.0316820.g005 In Fig 5, the events in the fault tree model are transformed into the nodes of the Bayesian network model, and the logical symbols are directly transformed into directed edges. Meanwhile, the dependency relationship between nodes is represented by a conditional probability table. Taking the "or" gate in the fault tree model as an example, its conditional probability table is shown in Table 2, and 0 and 1 respectively indicate the occurrence and non-occurrence of the events. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. The conditional probability distribution corresponding to the "or" gate. https://doi.org/10.1371/journal.pone.0316820.t002 The Bayesian network model uses probability integration to represent uncertainty, and the fundamental rule of probability integration is the probability of joint events, as shown in Formula (1). (1) Where P(B) is the occurrence probability of the event B; P (A, B) is the joint probability of the event A and B; P(A|B) is the conditional probability of the event A occurring under the condition that the event B occurs. By substituting the prior probability into the joint probability distribution, the occurrence probability of the highest-level event can be obtained, as shown in Formula (2). (2) Where Pa(Xi) is the occurrence probability of the disaster-causing factor Xi, which corresponds to the prior probability; P(Xi|Pa(Xi)) is the conditional probability of the higher-level factor occurring under the condition that the disaster-causing factor Xi occurs; P(D) is the occurrence probability of the highest-level event, which corresponds to the hazard probability of disaster events of gas pipelines in karst areas. 2.1.3. Determination of prior probability. At present, the failure database of gas pipelines in karst areas has not been established, and it is often difficult to ensure the accuracy and applicability to evaluate the occurrence probability of basic events by combining statistical analysis model. Therefore, it is necessary to determine the failure possibility of basic events based on expert experience. However, it is worth noting that even industry experts can only give some general judgments on the probability of fuzzy events that cause natural gas pipeline disasters (such as the lack of safety supervision, inadequate emergency support, etc.) [26]. In this regard, fuzzy mathematics is a good method for solving the problem of insufficient known data. In the theory of fuzzy mathematics, the index score in the evaluation system can reflect the occurrence possibility of an event. Assuming that the V0 is the maximum score, the score values of each evaluation index can be divided into five intervals, corresponding to different possibility levels [27], as shown in Table 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. The correspondence between score intervals and possibility levels. https://doi.org/10.1371/journal.pone.0316820.t003 Considering that the differences of experts’ knowledge level, experience level, information source and justice degree will lead to the differences of their decision-making, the analytic hierarchy process (AHP) model of expert authority evaluation is established, as shown in Fig 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. AHP model of expert authority evaluation. https://doi.org/10.1371/journal.pone.0316820.g006 Next, based on the judgment matrix, the influence weights of different related factors on the overall evaluation ability of experts, as well as the weights of different experts on the same evaluation ability related factors are investigated respectively. The construction principle of the evaluation matrix is that the cij in the matrix indicates the superiority of the indicator ui compared with the indicator uj, which conforms to the setting of Formula (3). (3) By revising the evaluation opinions based on the weights of experts’ comprehensive ability, it is helpful to make the evaluation results more in line with the objective reality. In this study, an improved AHP model based on the 1–9 exponential scale is used to determine the expert authority weights, and its specific scale value is shown in Table 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Exponential scale value. https://doi.org/10.1371/journal.pone.0316820.t004 Next, the weight vectors are used to represent the influence weights of different experts on the same evaluation ability related factor. Assuming that there are n experts in total, the cst represents the superiority of the s-th expert over the t-th expert for the i-th evaluation ability related factor, and its normalized value is c’st. In this regard, the weight Wis of the s-th expert relative to the i-th evaluation ability related factor is shown in Formula (4). (4) The basic principle of the AHP method is to compare the importance of influencing factors. Then, based on consistency testing, the evaluation results are investigated to determine if their inconsistency exceeds the allowable range. The fuzzy set theory uses fuzzy numbers to represent the relationship between uncertain variables and membership functions. In this study, the trapezoidal fuzzy number is selected to capture the uncertainty of experts’ judgment. In this regard, the risk levels of basic events corresponding to different indexs can be further transformed into the trapezoidal fuzzy numbers [28], as shown in Fig 7. For a certain fuzzy number M, its corresponding membership function model is shown in Formula (5). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Membership functions corresponding to trapezoidal the fuzzy numbers. https://doi.org/10.1371/journal.pone.0316820.g007 (5) Where fmax and fmin are the upper and lower limits of the fuzzy numbers, respectively. The α-cut set of fuzzy set is selected to comprehensively process the experts’ evaluation opinions, and then the corresponding fuzzy number function can be obtained. The fuzzy number function based on experts’ evaluation opinions is a fuzzy set within [0, 1]. It is necessary to convert it into a clear value, that is, fuzzy possibility score. It represents experts’ trust in the occurring of an basic event. Based on the left and right fuzzy sorting method, the fuzzy number can be converted into fuzzy possibility score [29]. The corresponding maximum fuzzy set S(f)max and minimum fuzzy set S(f)min are shown in Formula (6) and Formula (7) respectively. (6)(7) On this basis, the left and right fuzzy possibility score of the fuzzy number M are shown in Formula (8) and Formula (9) respectively. (8)(9) Where sup is the supremum of the set; EN(f) is the membership function corresponding to the fuzzy number M; ∧ is the arithmetic operation of taking the smaller of two values in fuzzy mathematics theory. On this basis, the fuzzy possibility score S of the fuzzy number M is shown in Formula (10). (10) Then, the fuzzy possibility score of each evaluation index is transformed into fuzzy failure probability. It is taken as the evaluation value of the prior probability of the corresponding disaster-causing factor, and its specific calculation is shown in Formulas (11) and (12). (11)(12) 2.2. Vulnerability analysis of gas pipeline The vulnerability refers to the possibility that the disaster-bearing body is susceptible to the damage or destruction caused by disaster-causing factors [30]. That is, the structure failure probability under different levels of disasters. In engineering practice, the vulnerability has been given different definitions and can be generally classified into the following three categories: (1) The possibility of potential destructive events causing damage to the disaster-bearing bodies; (2) The ability of the disaster-bearing bodies to predict, respond to, and handle disasters, as well as the ability to recover after disasters; (3) The comprehensive measurement of economic and social capacity of disaster-bearing bodies to respond to disaster events and disaster risks. The vulnerability level of gas pipelines is closely related to the disaster type and intensity, pipeline structure and function, the temporal and spatial configuration relationship between the disasters and the pipelines, and the ability of human society to prevent and reduce disasters [31]. 2.2.1. Construction of vulnerable function based on structural reliability theory. In the field of pipeline structural reliability research, the load applied to the pipelines that causes mechanical response is called generalized load, and the ability of the pipeline structure to resist generalized load is called generalized resistance [32]. The critical state function for pipeline structural damage is defined as D, which is expressed as the difference between the generalized resistance and the generalized load. When the function value is greater than 0, the pipeline structure is in a safe state; When the value is equal to 0, it is in a critical state; When the value is less than 0, it will enter the invalid state. In this regard, the basic principle equation of the structure failure probability Pf of gas pipelines is shown in Formula (13). (13) For the gas pipelines in karst areas, the vulnerability depends on the combined effects of the damage intensity of disaster load to the pipelines, and the ability of pipelines to resist disaster damage. The pipeline vulnerability probability corresponds to the structure failure probability in the structural reliability theory. Therefore, the Formula (13) is transformed to obtain the pipeline vulnerability function, as shown in Formula (14). (14) Where V is the vulnerability; Ir is the ability of pipelines to resist disaster damage; Id is the damage intensity of disaster load to the pipelines. 2.2.2. Selection of evaluation index. Karst collapse is a kind of karst dynamic geological effect and phenomenon, and its specific development characteristics are as follows: under the influences of natural and human activities, caves and cracks develop at the top of bedrock. Their continuous upward expansion will cause the deformation of the upper rock and soil, which eventually leads to the rupture of the karst roof and the formation of surface collapse pits. In order to ensure the stability of gas pipelines and reduce the impact of external factors, gas pipelines are usually laid underground. In this regard, the deformation of buried gas pipelines will be constrained by the surrounding soil, as shown in Fig 8. Meanwhile, the mechanical parameters and deformation performance of buried pipelines are significantly different from the surrounding soil, which will lead to complex relative deformation between the pipelines and the soil in the process of karst collapse. On this basis, the variables of vulnerability function are selected from three aspects: karst collapse parameters, pipeline operation parameters and soil properties parameters, as shown in Table 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Schematic diagram of pipeline deformation under the influence of karst collapse. https://doi.org/10.1371/journal.pone.0316820.g008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 5. The variables of vulnerability function. https://doi.org/10.1371/journal.pone.0316820.t005 Based on the engineering practice and previous research [33], the following influence laws are summarized: firstly, for the karst collapse parameters, the greater the length and the width of karst soil collapse, the greater the pipeline displacement caused by soil dislocation, and the greater the soil load on the pipelines. Under the action of external force, buried pipelines will produce certain plastic deformation, especially in internal and external defects and weak structures. The stability of buried pipelines is greatly affected, and the service life of pipelines will be significantly reduced if they continue to be used for a long time. The greater the thickness of karst overburden, the more stable the upper soil structure. In this regard, the influence of underground karst collapse on the buried pipelines will weaken in the process of upward conduction, so the mechanical response of pipelines under the impact of soil load will be reduced. Secondly, for the pipeline operation parameters, the greater the wall thickness, the smaller the diameter-thickness ratio, and the greater the pipeline stiffness. The increase of pipe buried depth is beneficial to strengthen the soil arching effect and reduce the karst collapse range, thus reducing its damage capacity. However, karst soil has a strong compressive effect on deeply buried gas pipelines, leading to a sharp increase in compressive stress. Meanwhile, deep buried gas pipelines are more easily affected by the dislocation of underground soil. With the increase of service time year by year, the aging phenomenon will lead to a significant decline in pipeline strength, so service time is a key parameter affecting the vulnerability of pipelines. Thirdly, for the soil property parameters, the increase of elastic modulus of soil will reduce the deformation ability of soil and increase the soil reaction force. Under the constraint of large soil reaction forces, the compression effect on the pipelines will increase, making it difficult to utilize the adaptability of gas pipeline ductility to soil dislocation. A small friction coefficient is beneficial for reducing the soil constraints on buried pipelines under the influence of karst collapse and minimizing mutual friction. In this regard, it can educe the stress concentration of dangerous sections of buried pipelines in karst collapse areas. For the selected eight variables, it is assumed that the critical parameters that lead to the failure of gas pipelines in karst areas are Kscl, c, Kscw, c, Kot, c, Pwt, c, Pbd, c, Sem, c, Fc, c and Lps, crespectively in a specific research scenario. When the value of a variable is equal to its corresponding critical parameter, the gas pipeline will be in a critical failure state. On this basis, the vulnerability parameters of different variables are defined as a1, a2, a3, a4, a5, a6, a7 and a8 respectively. Taking the length of karst soil collapse Kscl, c as an example, its vulnerability parameter is shown in Formula (15). Then, the vulnerability probability function can be further transformed into Formula (16). (15)(16) 2.2.3. Calculation of vulnerability probability. In this study, the vulnerability probability is calculated based on reliability evaluation method. The commonly used reliability calculation methods include second moment method, Monte Carlo method, etc. On the premise that the distribution of random variables is not clear, the first order second moment method uses the mean and standard deviation of random variables to perform linear Taylor series expansion on the functional function. The first order second moment method mainly includes the central point method and the design checkpoint method. The central point method is simple to calculate, but it is difficult to obtain the consistent reliability index for the functional functions with the same meaning but different mathematical expressions. In this regard, the design checkpoint point method can compensate for the deficiency. Meanwhile, the function expansion of the design checkpoint method at the average value is reasonable, so it can be chosen to calculate the vulnerability probability of gas pipelines in karst areas. Firstly, the reliability index β is introduced. Then, it is assumed that the vulnerability function V obeys a normal distribution, with an average value of μv and a standard deviation of σv. In this regard, the vulnerability probability function is shown in Formula (17). (17) If v = μv + σvu, then dv = σvdu. When v = 0, then u = -μv/σv. When v → -∞, then u → -∞. In this regard, the vulnerability probability function can be further transformed into Formula (18). (18) Where Φ(β) is the cumulative function of standard normal distribution; β is the pipeline reliability index, and β = μv/σv. For the general form of vulnerability function: (19) Where v1, v2, …, vn are n independent normal random variables, with an average of μj (j = 1, 2, …, n) and a standard deviation of σj (j = 1, 2, …, n). Assuming that a design checkpoint is represented as v* = (v1*, v2*, …, vn*), and g(v*) is equal to 0. Then, the functional function is expanded according to the Taylor series, and its first order expansion is shown in Formula (20). (20) The mean and standard deviation of ML are shown in Formula (21) and Formula (22). (21)(22) On this basis, the expression of reliability index β is shown in Formula (23). (23) The disaster-causing possibility of gas pipelines in karst areas is comprehensively investigated from two dimensions: the hazard of disaster events and the vulnerability of gas pipelines. If the disaster-causing probability is P, then P = P(D)P(V). According to the standard "API 581–2008 Risk-based Inspection Technology" issued by the American Petroleum Institute (API) and "DNV-RP-F116 Integrity Management of Submarine Pipeline Systems" issued by the Det Norske Veritas (DNV), the disaster-causing probability of gas pipelines in karst areas is divided into five grades, as shown in Table 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 6. Classification criteria for disaster-causing probability. https://doi.org/10.1371/journal.pone.0316820.t006 2.1. Hazard analysis of disaster events The hazard probability of disaster events is defined as the occurrence probability of all types of disaster events that threaten the safety of human beings and material wealth in the area. That is, the probability distribution of multiple risk inducing factors turning into serious disaster events [24]. 2.1.1. Establishment of hazard evaluation index system. Based on the disaster system theory, the hazard evaluation index system is established from three aspects: the activity of disaster-prone environment, the risk of disaster factors and the vulnerability of disaster-bearing bodies. (1) The disaster-prone environment. In a broad sense, the disaster-prone environment refers to a comprehensive social-natural system composed of multiple factors of atmosphere, lithosphere, biosphere, hydrosphere, and social material and cultural sphere. However, the formation of the disaster-prone environment is not simply the direct superposition of the related factors, but depends on the material circulation, energy flow and information transmission among various factors. The disaster-prone environment of gas pipeline disasters in karst areas is the breeding place of disaster factors, and it is also the indirect medium of the interaction between disaster factors and disaster-bearing bodies. The activity of disaster-prone environment (a) is closely related to the personnel unsafe state (k1), pipeline unsafe state (k2), environmental unsafe state (k3) and management loopholes (k4), as shown in Fig 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. The structure of disaster-prone environment. https://doi.org/10.1371/journal.pone.0316820.g001 Specifically, the personnel unsafe state (k1) mainly includes the low level of education (t1), the lack of professional skills (t2), the weak awareness of pipeline protection (t3), weak safety awareness (t4) and insufficient emergency response capacity (t5). The pipeline unsafe state (k2) mainly includes the pipeline aging (t6), pipeline corrosion (t7), the failure of safety facilities (t8) and construction technology defect (t9). The environmental unsafe state (k3) mainly includes the intensity of human activity (t10), construction and vibration (t11), groundwater extraction (t12), economic development level (t13), legal environment (t14), karst geological development (t15), groundwater activity (t16), overburden characteristics (t17), structural condition (t18), topographic features (t19) and meteorological condition (t20). The management loopholes (k4) mainly includes the lack of safety supervision (t21), the lack of safety publicity (t22), inadequate emergency support (t23) and unreasonable rules and regulations (t24). (2) The disaster factors. The disaster factors are the driving force and direct cause of hazards, which are manifested as various sudden extreme events in natural or social environments that threaten human life, cause material property losses, and cause ecological damage. They include both natural anomalies such as mudslides, tsunamis, tornadoes, and earthquakes, as well as anthropogenic anomalies such as desertification, hazardous spills, traffic accidents, and explosions. The disaster factors of gas pipelines in karst areas are the direct causes of casualties, economic losses and environmental damage, mainly including fire thermal radiation and explosion shock wave overpressure. The risk of disaster factors is mainly influenced by the probability, intensity, scope and duration of fire and explosion hazards. Among the various disaster factors caused by gas pipeline failure and leakage, the occurrence of a single disaster can lead to the subsequent occurrence of other secondary disasters. Therefore, the disaster factors of gas pipelines in karst areas show a chain formation process. The risk of disaster factors (b) is closely related to the fire risk (k5) and explosion risk (k6), as shown in Fig 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. The structure of disaster factors. https://doi.org/10.1371/journal.pone.0316820.g002 Specifically, the fire risk (k6) mainly includes the source of fire (t25), thermal radiation flux (t26), smoke and dust (t27), noxious fumes (t28), the types of gas leakage (t29), the scope of fire impact (t30), radiation time (t31), distance (t32) and space-limited situation (t36). The explosion risk (k5) mainly includes the shock wave overpressure (t33), noise (t34), the source of fire (t25), the types of gas leakage (t29), distance (t32), the scope of explosion impact (t35) and space-limited situation (t36). (3) The disaster-bearing body. The disaster-bearing bodies refers to the direct natural and social subjects affected by disasters, including various social resources such as human life, material property, social production and operation systems, as well as various natural resources such as animals, plants, soil, water sources, and atmosphere. The disaster-bearing bodies of gas pipeline disasters in karst areas are the direct object that bears the hazard of disaster factors, which mainly includes social, economic and natural disaster-bearing bodies along the gas pipelines. Different disaster-bearing bodies are closely related and interact with each other. The vulnerability of disaster-bearing bodies is mainly used to measure the difficulty of their being damaged by disasters. On the one hand, it depends on the physical characteristics of the disaster-bearing bodies, such as the strength and ductility level of gas pipelines. On the other hand, it depends on the cultivated disaster-bearing ability, that is, the multiple abilities of resistance, rescue and recovery in the face of disasters. The vulnerability of disaster-bearing bodies (c) is closely related to the social disaster-bearing body (k7), economic disaster-bearing body (k8), natural disaster-bearing body (k9) and disaster bearing capacity (k10), as shown in Fig 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. The structure of disaster-bearing bodies. https://doi.org/10.1371/journal.pone.0316820.g003 Specifically, the social disaster-bearing body (k7) mainly includes the human casualties (t37), the destruction of buildings and structures (t38) and the destruction of lifeline engineering (t39). The economic disaster-bearing body (k8) mainly includes the material property loss (t40), economic compensation (t41), resource waste (t42) and production and business suspension (t43). The natural disaster-bearing body (k9) mainly includes the atmospheric pollution (t44), soil hardening (t45), the damage to vegetation (t46) and the harm to animals (t47). The disaster bearing capacity (k10) mainly includes the early warning capability (t48), disaster resistance capability (t49), disaster relief capacity (t50) and recovery capacity (t51). 2.1.2. Establishment of Bayesian network model. The Bayesian network is a directed acyclic graph based on graph theory and probability theory, which is mainly used to describe the dependency relationships and probability distributions of a group of random variables. It consists of nodes, directed edges and probabilities, and a simple Bayesian network structure is shown in Fig 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. A simple Bayesian network structure. https://doi.org/10.1371/journal.pone.0316820.g004 In Fig 4, nodes A, B1, B2, C1 and C2 represent a group of random variables, and the arrows indicate their dependency relationships. The C1 points to B1, indicating that the C1 is the parent node of the B1, and the B1 is the child node of the C1. Meanwhile, the C1 and C2 have no parent node, so they are root nodes. The root nodes have corresponding prior probability distribution, and non-root nodes have corresponding conditional probability distribution. The Bayesian network can update the failure probability of the system by adjusting the probability distribution of nodes. Taking Fig 4 as an example, the adjustment process of probability distribution of nodes in the Bayesian network is investigated: assuming that at a certain moment t, the probability that the occurrence of C1 causes the occurrence of B1 is 0.2, the probability that the occurrence of C2 causes the occurrence of B1 is 0.3. The probability of B1 occurring when C1 and C2 both occur is 0.5, and the probability of B1 occurring when neither C1 and C2 occurs is 0.1. In this regard, the conditional probability distribution corresponding to the node B1 is shown in Table 1, and 0 and 1 respectively indicate the occurrence and non-occurrence of the root nodes. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. The conditional probability distribution corresponding to the node B1. https://doi.org/10.1371/journal.pone.0316820.t001 In the Bayesian network model, the probability calculation of event occurrence can be equivalent to the propagation and updating of beliefs [25]. In this regard, the calculation process is simple and the results are accurate. Therefore, the Bayesian network model can be used to evaluate the hazard probability of disaster events of gas pipelines in karst areas. Meanwhile, the fault tree model and Bayesian network model have similarities in principle, and the deductive reasoning function of the fault tree model can quickly determine the dependency relationships between multiple events. In this regard, the Bayesian network model can be established based on the corresponding mapping relationship. The transformation of the “or” gate and “and” gate in the fault tree model to the Bayesian network model is shown in Fig 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Schematic diagram of transformation between the fault tree model and the Bayesian network model. The “or” gate (A), and“and” gate (B). https://doi.org/10.1371/journal.pone.0316820.g005 In Fig 5, the events in the fault tree model are transformed into the nodes of the Bayesian network model, and the logical symbols are directly transformed into directed edges. Meanwhile, the dependency relationship between nodes is represented by a conditional probability table. Taking the "or" gate in the fault tree model as an example, its conditional probability table is shown in Table 2, and 0 and 1 respectively indicate the occurrence and non-occurrence of the events. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. The conditional probability distribution corresponding to the "or" gate. https://doi.org/10.1371/journal.pone.0316820.t002 The Bayesian network model uses probability integration to represent uncertainty, and the fundamental rule of probability integration is the probability of joint events, as shown in Formula (1). (1) Where P(B) is the occurrence probability of the event B; P (A, B) is the joint probability of the event A and B; P(A|B) is the conditional probability of the event A occurring under the condition that the event B occurs. By substituting the prior probability into the joint probability distribution, the occurrence probability of the highest-level event can be obtained, as shown in Formula (2). (2) Where Pa(Xi) is the occurrence probability of the disaster-causing factor Xi, which corresponds to the prior probability; P(Xi|Pa(Xi)) is the conditional probability of the higher-level factor occurring under the condition that the disaster-causing factor Xi occurs; P(D) is the occurrence probability of the highest-level event, which corresponds to the hazard probability of disaster events of gas pipelines in karst areas. 2.1.3. Determination of prior probability. At present, the failure database of gas pipelines in karst areas has not been established, and it is often difficult to ensure the accuracy and applicability to evaluate the occurrence probability of basic events by combining statistical analysis model. Therefore, it is necessary to determine the failure possibility of basic events based on expert experience. However, it is worth noting that even industry experts can only give some general judgments on the probability of fuzzy events that cause natural gas pipeline disasters (such as the lack of safety supervision, inadequate emergency support, etc.) [26]. In this regard, fuzzy mathematics is a good method for solving the problem of insufficient known data. In the theory of fuzzy mathematics, the index score in the evaluation system can reflect the occurrence possibility of an event. Assuming that the V0 is the maximum score, the score values of each evaluation index can be divided into five intervals, corresponding to different possibility levels [27], as shown in Table 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. The correspondence between score intervals and possibility levels. https://doi.org/10.1371/journal.pone.0316820.t003 Considering that the differences of experts’ knowledge level, experience level, information source and justice degree will lead to the differences of their decision-making, the analytic hierarchy process (AHP) model of expert authority evaluation is established, as shown in Fig 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. AHP model of expert authority evaluation. https://doi.org/10.1371/journal.pone.0316820.g006 Next, based on the judgment matrix, the influence weights of different related factors on the overall evaluation ability of experts, as well as the weights of different experts on the same evaluation ability related factors are investigated respectively. The construction principle of the evaluation matrix is that the cij in the matrix indicates the superiority of the indicator ui compared with the indicator uj, which conforms to the setting of Formula (3). (3) By revising the evaluation opinions based on the weights of experts’ comprehensive ability, it is helpful to make the evaluation results more in line with the objective reality. In this study, an improved AHP model based on the 1–9 exponential scale is used to determine the expert authority weights, and its specific scale value is shown in Table 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Exponential scale value. https://doi.org/10.1371/journal.pone.0316820.t004 Next, the weight vectors are used to represent the influence weights of different experts on the same evaluation ability related factor. Assuming that there are n experts in total, the cst represents the superiority of the s-th expert over the t-th expert for the i-th evaluation ability related factor, and its normalized value is c’st. In this regard, the weight Wis of the s-th expert relative to the i-th evaluation ability related factor is shown in Formula (4). (4) The basic principle of the AHP method is to compare the importance of influencing factors. Then, based on consistency testing, the evaluation results are investigated to determine if their inconsistency exceeds the allowable range. The fuzzy set theory uses fuzzy numbers to represent the relationship between uncertain variables and membership functions. In this study, the trapezoidal fuzzy number is selected to capture the uncertainty of experts’ judgment. In this regard, the risk levels of basic events corresponding to different indexs can be further transformed into the trapezoidal fuzzy numbers [28], as shown in Fig 7. For a certain fuzzy number M, its corresponding membership function model is shown in Formula (5). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Membership functions corresponding to trapezoidal the fuzzy numbers. https://doi.org/10.1371/journal.pone.0316820.g007 (5) Where fmax and fmin are the upper and lower limits of the fuzzy numbers, respectively. The α-cut set of fuzzy set is selected to comprehensively process the experts’ evaluation opinions, and then the corresponding fuzzy number function can be obtained. The fuzzy number function based on experts’ evaluation opinions is a fuzzy set within [0, 1]. It is necessary to convert it into a clear value, that is, fuzzy possibility score. It represents experts’ trust in the occurring of an basic event. Based on the left and right fuzzy sorting method, the fuzzy number can be converted into fuzzy possibility score [29]. The corresponding maximum fuzzy set S(f)max and minimum fuzzy set S(f)min are shown in Formula (6) and Formula (7) respectively. (6)(7) On this basis, the left and right fuzzy possibility score of the fuzzy number M are shown in Formula (8) and Formula (9) respectively. (8)(9) Where sup is the supremum of the set; EN(f) is the membership function corresponding to the fuzzy number M; ∧ is the arithmetic operation of taking the smaller of two values in fuzzy mathematics theory. On this basis, the fuzzy possibility score S of the fuzzy number M is shown in Formula (10). (10) Then, the fuzzy possibility score of each evaluation index is transformed into fuzzy failure probability. It is taken as the evaluation value of the prior probability of the corresponding disaster-causing factor, and its specific calculation is shown in Formulas (11) and (12). (11)(12) 2.1.1. Establishment of hazard evaluation index system. Based on the disaster system theory, the hazard evaluation index system is established from three aspects: the activity of disaster-prone environment, the risk of disaster factors and the vulnerability of disaster-bearing bodies. (1) The disaster-prone environment. In a broad sense, the disaster-prone environment refers to a comprehensive social-natural system composed of multiple factors of atmosphere, lithosphere, biosphere, hydrosphere, and social material and cultural sphere. However, the formation of the disaster-prone environment is not simply the direct superposition of the related factors, but depends on the material circulation, energy flow and information transmission among various factors. The disaster-prone environment of gas pipeline disasters in karst areas is the breeding place of disaster factors, and it is also the indirect medium of the interaction between disaster factors and disaster-bearing bodies. The activity of disaster-prone environment (a) is closely related to the personnel unsafe state (k1), pipeline unsafe state (k2), environmental unsafe state (k3) and management loopholes (k4), as shown in Fig 1. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 1. The structure of disaster-prone environment. https://doi.org/10.1371/journal.pone.0316820.g001 Specifically, the personnel unsafe state (k1) mainly includes the low level of education (t1), the lack of professional skills (t2), the weak awareness of pipeline protection (t3), weak safety awareness (t4) and insufficient emergency response capacity (t5). The pipeline unsafe state (k2) mainly includes the pipeline aging (t6), pipeline corrosion (t7), the failure of safety facilities (t8) and construction technology defect (t9). The environmental unsafe state (k3) mainly includes the intensity of human activity (t10), construction and vibration (t11), groundwater extraction (t12), economic development level (t13), legal environment (t14), karst geological development (t15), groundwater activity (t16), overburden characteristics (t17), structural condition (t18), topographic features (t19) and meteorological condition (t20). The management loopholes (k4) mainly includes the lack of safety supervision (t21), the lack of safety publicity (t22), inadequate emergency support (t23) and unreasonable rules and regulations (t24). (2) The disaster factors. The disaster factors are the driving force and direct cause of hazards, which are manifested as various sudden extreme events in natural or social environments that threaten human life, cause material property losses, and cause ecological damage. They include both natural anomalies such as mudslides, tsunamis, tornadoes, and earthquakes, as well as anthropogenic anomalies such as desertification, hazardous spills, traffic accidents, and explosions. The disaster factors of gas pipelines in karst areas are the direct causes of casualties, economic losses and environmental damage, mainly including fire thermal radiation and explosion shock wave overpressure. The risk of disaster factors is mainly influenced by the probability, intensity, scope and duration of fire and explosion hazards. Among the various disaster factors caused by gas pipeline failure and leakage, the occurrence of a single disaster can lead to the subsequent occurrence of other secondary disasters. Therefore, the disaster factors of gas pipelines in karst areas show a chain formation process. The risk of disaster factors (b) is closely related to the fire risk (k5) and explosion risk (k6), as shown in Fig 2. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 2. The structure of disaster factors. https://doi.org/10.1371/journal.pone.0316820.g002 Specifically, the fire risk (k6) mainly includes the source of fire (t25), thermal radiation flux (t26), smoke and dust (t27), noxious fumes (t28), the types of gas leakage (t29), the scope of fire impact (t30), radiation time (t31), distance (t32) and space-limited situation (t36). The explosion risk (k5) mainly includes the shock wave overpressure (t33), noise (t34), the source of fire (t25), the types of gas leakage (t29), distance (t32), the scope of explosion impact (t35) and space-limited situation (t36). (3) The disaster-bearing body. The disaster-bearing bodies refers to the direct natural and social subjects affected by disasters, including various social resources such as human life, material property, social production and operation systems, as well as various natural resources such as animals, plants, soil, water sources, and atmosphere. The disaster-bearing bodies of gas pipeline disasters in karst areas are the direct object that bears the hazard of disaster factors, which mainly includes social, economic and natural disaster-bearing bodies along the gas pipelines. Different disaster-bearing bodies are closely related and interact with each other. The vulnerability of disaster-bearing bodies is mainly used to measure the difficulty of their being damaged by disasters. On the one hand, it depends on the physical characteristics of the disaster-bearing bodies, such as the strength and ductility level of gas pipelines. On the other hand, it depends on the cultivated disaster-bearing ability, that is, the multiple abilities of resistance, rescue and recovery in the face of disasters. The vulnerability of disaster-bearing bodies (c) is closely related to the social disaster-bearing body (k7), economic disaster-bearing body (k8), natural disaster-bearing body (k9) and disaster bearing capacity (k10), as shown in Fig 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 3. The structure of disaster-bearing bodies. https://doi.org/10.1371/journal.pone.0316820.g003 Specifically, the social disaster-bearing body (k7) mainly includes the human casualties (t37), the destruction of buildings and structures (t38) and the destruction of lifeline engineering (t39). The economic disaster-bearing body (k8) mainly includes the material property loss (t40), economic compensation (t41), resource waste (t42) and production and business suspension (t43). The natural disaster-bearing body (k9) mainly includes the atmospheric pollution (t44), soil hardening (t45), the damage to vegetation (t46) and the harm to animals (t47). The disaster bearing capacity (k10) mainly includes the early warning capability (t48), disaster resistance capability (t49), disaster relief capacity (t50) and recovery capacity (t51). 2.1.2. Establishment of Bayesian network model. The Bayesian network is a directed acyclic graph based on graph theory and probability theory, which is mainly used to describe the dependency relationships and probability distributions of a group of random variables. It consists of nodes, directed edges and probabilities, and a simple Bayesian network structure is shown in Fig 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 4. A simple Bayesian network structure. https://doi.org/10.1371/journal.pone.0316820.g004 In Fig 4, nodes A, B1, B2, C1 and C2 represent a group of random variables, and the arrows indicate their dependency relationships. The C1 points to B1, indicating that the C1 is the parent node of the B1, and the B1 is the child node of the C1. Meanwhile, the C1 and C2 have no parent node, so they are root nodes. The root nodes have corresponding prior probability distribution, and non-root nodes have corresponding conditional probability distribution. The Bayesian network can update the failure probability of the system by adjusting the probability distribution of nodes. Taking Fig 4 as an example, the adjustment process of probability distribution of nodes in the Bayesian network is investigated: assuming that at a certain moment t, the probability that the occurrence of C1 causes the occurrence of B1 is 0.2, the probability that the occurrence of C2 causes the occurrence of B1 is 0.3. The probability of B1 occurring when C1 and C2 both occur is 0.5, and the probability of B1 occurring when neither C1 and C2 occurs is 0.1. In this regard, the conditional probability distribution corresponding to the node B1 is shown in Table 1, and 0 and 1 respectively indicate the occurrence and non-occurrence of the root nodes. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 1. The conditional probability distribution corresponding to the node B1. https://doi.org/10.1371/journal.pone.0316820.t001 In the Bayesian network model, the probability calculation of event occurrence can be equivalent to the propagation and updating of beliefs [25]. In this regard, the calculation process is simple and the results are accurate. Therefore, the Bayesian network model can be used to evaluate the hazard probability of disaster events of gas pipelines in karst areas. Meanwhile, the fault tree model and Bayesian network model have similarities in principle, and the deductive reasoning function of the fault tree model can quickly determine the dependency relationships between multiple events. In this regard, the Bayesian network model can be established based on the corresponding mapping relationship. The transformation of the “or” gate and “and” gate in the fault tree model to the Bayesian network model is shown in Fig 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 5. Schematic diagram of transformation between the fault tree model and the Bayesian network model. The “or” gate (A), and“and” gate (B). https://doi.org/10.1371/journal.pone.0316820.g005 In Fig 5, the events in the fault tree model are transformed into the nodes of the Bayesian network model, and the logical symbols are directly transformed into directed edges. Meanwhile, the dependency relationship between nodes is represented by a conditional probability table. Taking the "or" gate in the fault tree model as an example, its conditional probability table is shown in Table 2, and 0 and 1 respectively indicate the occurrence and non-occurrence of the events. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 2. The conditional probability distribution corresponding to the "or" gate. https://doi.org/10.1371/journal.pone.0316820.t002 The Bayesian network model uses probability integration to represent uncertainty, and the fundamental rule of probability integration is the probability of joint events, as shown in Formula (1). (1) Where P(B) is the occurrence probability of the event B; P (A, B) is the joint probability of the event A and B; P(A|B) is the conditional probability of the event A occurring under the condition that the event B occurs. By substituting the prior probability into the joint probability distribution, the occurrence probability of the highest-level event can be obtained, as shown in Formula (2). (2) Where Pa(Xi) is the occurrence probability of the disaster-causing factor Xi, which corresponds to the prior probability; P(Xi|Pa(Xi)) is the conditional probability of the higher-level factor occurring under the condition that the disaster-causing factor Xi occurs; P(D) is the occurrence probability of the highest-level event, which corresponds to the hazard probability of disaster events of gas pipelines in karst areas. 2.1.3. Determination of prior probability. At present, the failure database of gas pipelines in karst areas has not been established, and it is often difficult to ensure the accuracy and applicability to evaluate the occurrence probability of basic events by combining statistical analysis model. Therefore, it is necessary to determine the failure possibility of basic events based on expert experience. However, it is worth noting that even industry experts can only give some general judgments on the probability of fuzzy events that cause natural gas pipeline disasters (such as the lack of safety supervision, inadequate emergency support, etc.) [26]. In this regard, fuzzy mathematics is a good method for solving the problem of insufficient known data. In the theory of fuzzy mathematics, the index score in the evaluation system can reflect the occurrence possibility of an event. Assuming that the V0 is the maximum score, the score values of each evaluation index can be divided into five intervals, corresponding to different possibility levels [27], as shown in Table 3. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 3. The correspondence between score intervals and possibility levels. https://doi.org/10.1371/journal.pone.0316820.t003 Considering that the differences of experts’ knowledge level, experience level, information source and justice degree will lead to the differences of their decision-making, the analytic hierarchy process (AHP) model of expert authority evaluation is established, as shown in Fig 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 6. AHP model of expert authority evaluation. https://doi.org/10.1371/journal.pone.0316820.g006 Next, based on the judgment matrix, the influence weights of different related factors on the overall evaluation ability of experts, as well as the weights of different experts on the same evaluation ability related factors are investigated respectively. The construction principle of the evaluation matrix is that the cij in the matrix indicates the superiority of the indicator ui compared with the indicator uj, which conforms to the setting of Formula (3). (3) By revising the evaluation opinions based on the weights of experts’ comprehensive ability, it is helpful to make the evaluation results more in line with the objective reality. In this study, an improved AHP model based on the 1–9 exponential scale is used to determine the expert authority weights, and its specific scale value is shown in Table 4. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 4. Exponential scale value. https://doi.org/10.1371/journal.pone.0316820.t004 Next, the weight vectors are used to represent the influence weights of different experts on the same evaluation ability related factor. Assuming that there are n experts in total, the cst represents the superiority of the s-th expert over the t-th expert for the i-th evaluation ability related factor, and its normalized value is c’st. In this regard, the weight Wis of the s-th expert relative to the i-th evaluation ability related factor is shown in Formula (4). (4) The basic principle of the AHP method is to compare the importance of influencing factors. Then, based on consistency testing, the evaluation results are investigated to determine if their inconsistency exceeds the allowable range. The fuzzy set theory uses fuzzy numbers to represent the relationship between uncertain variables and membership functions. In this study, the trapezoidal fuzzy number is selected to capture the uncertainty of experts’ judgment. In this regard, the risk levels of basic events corresponding to different indexs can be further transformed into the trapezoidal fuzzy numbers [28], as shown in Fig 7. For a certain fuzzy number M, its corresponding membership function model is shown in Formula (5). Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 7. Membership functions corresponding to trapezoidal the fuzzy numbers. https://doi.org/10.1371/journal.pone.0316820.g007 (5) Where fmax and fmin are the upper and lower limits of the fuzzy numbers, respectively. The α-cut set of fuzzy set is selected to comprehensively process the experts’ evaluation opinions, and then the corresponding fuzzy number function can be obtained. The fuzzy number function based on experts’ evaluation opinions is a fuzzy set within [0, 1]. It is necessary to convert it into a clear value, that is, fuzzy possibility score. It represents experts’ trust in the occurring of an basic event. Based on the left and right fuzzy sorting method, the fuzzy number can be converted into fuzzy possibility score [29]. The corresponding maximum fuzzy set S(f)max and minimum fuzzy set S(f)min are shown in Formula (6) and Formula (7) respectively. (6)(7) On this basis, the left and right fuzzy possibility score of the fuzzy number M are shown in Formula (8) and Formula (9) respectively. (8)(9) Where sup is the supremum of the set; EN(f) is the membership function corresponding to the fuzzy number M; ∧ is the arithmetic operation of taking the smaller of two values in fuzzy mathematics theory. On this basis, the fuzzy possibility score S of the fuzzy number M is shown in Formula (10). (10) Then, the fuzzy possibility score of each evaluation index is transformed into fuzzy failure probability. It is taken as the evaluation value of the prior probability of the corresponding disaster-causing factor, and its specific calculation is shown in Formulas (11) and (12). (11)(12) 2.2. Vulnerability analysis of gas pipeline The vulnerability refers to the possibility that the disaster-bearing body is susceptible to the damage or destruction caused by disaster-causing factors [30]. That is, the structure failure probability under different levels of disasters. In engineering practice, the vulnerability has been given different definitions and can be generally classified into the following three categories: (1) The possibility of potential destructive events causing damage to the disaster-bearing bodies; (2) The ability of the disaster-bearing bodies to predict, respond to, and handle disasters, as well as the ability to recover after disasters; (3) The comprehensive measurement of economic and social capacity of disaster-bearing bodies to respond to disaster events and disaster risks. The vulnerability level of gas pipelines is closely related to the disaster type and intensity, pipeline structure and function, the temporal and spatial configuration relationship between the disasters and the pipelines, and the ability of human society to prevent and reduce disasters [31]. 2.2.1. Construction of vulnerable function based on structural reliability theory. In the field of pipeline structural reliability research, the load applied to the pipelines that causes mechanical response is called generalized load, and the ability of the pipeline structure to resist generalized load is called generalized resistance [32]. The critical state function for pipeline structural damage is defined as D, which is expressed as the difference between the generalized resistance and the generalized load. When the function value is greater than 0, the pipeline structure is in a safe state; When the value is equal to 0, it is in a critical state; When the value is less than 0, it will enter the invalid state. In this regard, the basic principle equation of the structure failure probability Pf of gas pipelines is shown in Formula (13). (13) For the gas pipelines in karst areas, the vulnerability depends on the combined effects of the damage intensity of disaster load to the pipelines, and the ability of pipelines to resist disaster damage. The pipeline vulnerability probability corresponds to the structure failure probability in the structural reliability theory. Therefore, the Formula (13) is transformed to obtain the pipeline vulnerability function, as shown in Formula (14). (14) Where V is the vulnerability; Ir is the ability of pipelines to resist disaster damage; Id is the damage intensity of disaster load to the pipelines. 2.2.2. Selection of evaluation index. Karst collapse is a kind of karst dynamic geological effect and phenomenon, and its specific development characteristics are as follows: under the influences of natural and human activities, caves and cracks develop at the top of bedrock. Their continuous upward expansion will cause the deformation of the upper rock and soil, which eventually leads to the rupture of the karst roof and the formation of surface collapse pits. In order to ensure the stability of gas pipelines and reduce the impact of external factors, gas pipelines are usually laid underground. In this regard, the deformation of buried gas pipelines will be constrained by the surrounding soil, as shown in Fig 8. Meanwhile, the mechanical parameters and deformation performance of buried pipelines are significantly different from the surrounding soil, which will lead to complex relative deformation between the pipelines and the soil in the process of karst collapse. On this basis, the variables of vulnerability function are selected from three aspects: karst collapse parameters, pipeline operation parameters and soil properties parameters, as shown in Table 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Schematic diagram of pipeline deformation under the influence of karst collapse. https://doi.org/10.1371/journal.pone.0316820.g008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 5. The variables of vulnerability function. https://doi.org/10.1371/journal.pone.0316820.t005 Based on the engineering practice and previous research [33], the following influence laws are summarized: firstly, for the karst collapse parameters, the greater the length and the width of karst soil collapse, the greater the pipeline displacement caused by soil dislocation, and the greater the soil load on the pipelines. Under the action of external force, buried pipelines will produce certain plastic deformation, especially in internal and external defects and weak structures. The stability of buried pipelines is greatly affected, and the service life of pipelines will be significantly reduced if they continue to be used for a long time. The greater the thickness of karst overburden, the more stable the upper soil structure. In this regard, the influence of underground karst collapse on the buried pipelines will weaken in the process of upward conduction, so the mechanical response of pipelines under the impact of soil load will be reduced. Secondly, for the pipeline operation parameters, the greater the wall thickness, the smaller the diameter-thickness ratio, and the greater the pipeline stiffness. The increase of pipe buried depth is beneficial to strengthen the soil arching effect and reduce the karst collapse range, thus reducing its damage capacity. However, karst soil has a strong compressive effect on deeply buried gas pipelines, leading to a sharp increase in compressive stress. Meanwhile, deep buried gas pipelines are more easily affected by the dislocation of underground soil. With the increase of service time year by year, the aging phenomenon will lead to a significant decline in pipeline strength, so service time is a key parameter affecting the vulnerability of pipelines. Thirdly, for the soil property parameters, the increase of elastic modulus of soil will reduce the deformation ability of soil and increase the soil reaction force. Under the constraint of large soil reaction forces, the compression effect on the pipelines will increase, making it difficult to utilize the adaptability of gas pipeline ductility to soil dislocation. A small friction coefficient is beneficial for reducing the soil constraints on buried pipelines under the influence of karst collapse and minimizing mutual friction. In this regard, it can educe the stress concentration of dangerous sections of buried pipelines in karst collapse areas. For the selected eight variables, it is assumed that the critical parameters that lead to the failure of gas pipelines in karst areas are Kscl, c, Kscw, c, Kot, c, Pwt, c, Pbd, c, Sem, c, Fc, c and Lps, crespectively in a specific research scenario. When the value of a variable is equal to its corresponding critical parameter, the gas pipeline will be in a critical failure state. On this basis, the vulnerability parameters of different variables are defined as a1, a2, a3, a4, a5, a6, a7 and a8 respectively. Taking the length of karst soil collapse Kscl, c as an example, its vulnerability parameter is shown in Formula (15). Then, the vulnerability probability function can be further transformed into Formula (16). (15)(16) 2.2.3. Calculation of vulnerability probability. In this study, the vulnerability probability is calculated based on reliability evaluation method. The commonly used reliability calculation methods include second moment method, Monte Carlo method, etc. On the premise that the distribution of random variables is not clear, the first order second moment method uses the mean and standard deviation of random variables to perform linear Taylor series expansion on the functional function. The first order second moment method mainly includes the central point method and the design checkpoint method. The central point method is simple to calculate, but it is difficult to obtain the consistent reliability index for the functional functions with the same meaning but different mathematical expressions. In this regard, the design checkpoint point method can compensate for the deficiency. Meanwhile, the function expansion of the design checkpoint method at the average value is reasonable, so it can be chosen to calculate the vulnerability probability of gas pipelines in karst areas. Firstly, the reliability index β is introduced. Then, it is assumed that the vulnerability function V obeys a normal distribution, with an average value of μv and a standard deviation of σv. In this regard, the vulnerability probability function is shown in Formula (17). (17) If v = μv + σvu, then dv = σvdu. When v = 0, then u = -μv/σv. When v → -∞, then u → -∞. In this regard, the vulnerability probability function can be further transformed into Formula (18). (18) Where Φ(β) is the cumulative function of standard normal distribution; β is the pipeline reliability index, and β = μv/σv. For the general form of vulnerability function: (19) Where v1, v2, …, vn are n independent normal random variables, with an average of μj (j = 1, 2, …, n) and a standard deviation of σj (j = 1, 2, …, n). Assuming that a design checkpoint is represented as v* = (v1*, v2*, …, vn*), and g(v*) is equal to 0. Then, the functional function is expanded according to the Taylor series, and its first order expansion is shown in Formula (20). (20) The mean and standard deviation of ML are shown in Formula (21) and Formula (22). (21)(22) On this basis, the expression of reliability index β is shown in Formula (23). (23) The disaster-causing possibility of gas pipelines in karst areas is comprehensively investigated from two dimensions: the hazard of disaster events and the vulnerability of gas pipelines. If the disaster-causing probability is P, then P = P(D)P(V). According to the standard "API 581–2008 Risk-based Inspection Technology" issued by the American Petroleum Institute (API) and "DNV-RP-F116 Integrity Management of Submarine Pipeline Systems" issued by the Det Norske Veritas (DNV), the disaster-causing probability of gas pipelines in karst areas is divided into five grades, as shown in Table 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 6. Classification criteria for disaster-causing probability. https://doi.org/10.1371/journal.pone.0316820.t006 2.2.1. Construction of vulnerable function based on structural reliability theory. In the field of pipeline structural reliability research, the load applied to the pipelines that causes mechanical response is called generalized load, and the ability of the pipeline structure to resist generalized load is called generalized resistance [32]. The critical state function for pipeline structural damage is defined as D, which is expressed as the difference between the generalized resistance and the generalized load. When the function value is greater than 0, the pipeline structure is in a safe state; When the value is equal to 0, it is in a critical state; When the value is less than 0, it will enter the invalid state. In this regard, the basic principle equation of the structure failure probability Pf of gas pipelines is shown in Formula (13). (13) For the gas pipelines in karst areas, the vulnerability depends on the combined effects of the damage intensity of disaster load to the pipelines, and the ability of pipelines to resist disaster damage. The pipeline vulnerability probability corresponds to the structure failure probability in the structural reliability theory. Therefore, the Formula (13) is transformed to obtain the pipeline vulnerability function, as shown in Formula (14). (14) Where V is the vulnerability; Ir is the ability of pipelines to resist disaster damage; Id is the damage intensity of disaster load to the pipelines. 2.2.2. Selection of evaluation index. Karst collapse is a kind of karst dynamic geological effect and phenomenon, and its specific development characteristics are as follows: under the influences of natural and human activities, caves and cracks develop at the top of bedrock. Their continuous upward expansion will cause the deformation of the upper rock and soil, which eventually leads to the rupture of the karst roof and the formation of surface collapse pits. In order to ensure the stability of gas pipelines and reduce the impact of external factors, gas pipelines are usually laid underground. In this regard, the deformation of buried gas pipelines will be constrained by the surrounding soil, as shown in Fig 8. Meanwhile, the mechanical parameters and deformation performance of buried pipelines are significantly different from the surrounding soil, which will lead to complex relative deformation between the pipelines and the soil in the process of karst collapse. On this basis, the variables of vulnerability function are selected from three aspects: karst collapse parameters, pipeline operation parameters and soil properties parameters, as shown in Table 5. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 8. Schematic diagram of pipeline deformation under the influence of karst collapse. https://doi.org/10.1371/journal.pone.0316820.g008 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 5. The variables of vulnerability function. https://doi.org/10.1371/journal.pone.0316820.t005 Based on the engineering practice and previous research [33], the following influence laws are summarized: firstly, for the karst collapse parameters, the greater the length and the width of karst soil collapse, the greater the pipeline displacement caused by soil dislocation, and the greater the soil load on the pipelines. Under the action of external force, buried pipelines will produce certain plastic deformation, especially in internal and external defects and weak structures. The stability of buried pipelines is greatly affected, and the service life of pipelines will be significantly reduced if they continue to be used for a long time. The greater the thickness of karst overburden, the more stable the upper soil structure. In this regard, the influence of underground karst collapse on the buried pipelines will weaken in the process of upward conduction, so the mechanical response of pipelines under the impact of soil load will be reduced. Secondly, for the pipeline operation parameters, the greater the wall thickness, the smaller the diameter-thickness ratio, and the greater the pipeline stiffness. The increase of pipe buried depth is beneficial to strengthen the soil arching effect and reduce the karst collapse range, thus reducing its damage capacity. However, karst soil has a strong compressive effect on deeply buried gas pipelines, leading to a sharp increase in compressive stress. Meanwhile, deep buried gas pipelines are more easily affected by the dislocation of underground soil. With the increase of service time year by year, the aging phenomenon will lead to a significant decline in pipeline strength, so service time is a key parameter affecting the vulnerability of pipelines. Thirdly, for the soil property parameters, the increase of elastic modulus of soil will reduce the deformation ability of soil and increase the soil reaction force. Under the constraint of large soil reaction forces, the compression effect on the pipelines will increase, making it difficult to utilize the adaptability of gas pipeline ductility to soil dislocation. A small friction coefficient is beneficial for reducing the soil constraints on buried pipelines under the influence of karst collapse and minimizing mutual friction. In this regard, it can educe the stress concentration of dangerous sections of buried pipelines in karst collapse areas. For the selected eight variables, it is assumed that the critical parameters that lead to the failure of gas pipelines in karst areas are Kscl, c, Kscw, c, Kot, c, Pwt, c, Pbd, c, Sem, c, Fc, c and Lps, crespectively in a specific research scenario. When the value of a variable is equal to its corresponding critical parameter, the gas pipeline will be in a critical failure state. On this basis, the vulnerability parameters of different variables are defined as a1, a2, a3, a4, a5, a6, a7 and a8 respectively. Taking the length of karst soil collapse Kscl, c as an example, its vulnerability parameter is shown in Formula (15). Then, the vulnerability probability function can be further transformed into Formula (16). (15)(16) 2.2.3. Calculation of vulnerability probability. In this study, the vulnerability probability is calculated based on reliability evaluation method. The commonly used reliability calculation methods include second moment method, Monte Carlo method, etc. On the premise that the distribution of random variables is not clear, the first order second moment method uses the mean and standard deviation of random variables to perform linear Taylor series expansion on the functional function. The first order second moment method mainly includes the central point method and the design checkpoint method. The central point method is simple to calculate, but it is difficult to obtain the consistent reliability index for the functional functions with the same meaning but different mathematical expressions. In this regard, the design checkpoint point method can compensate for the deficiency. Meanwhile, the function expansion of the design checkpoint method at the average value is reasonable, so it can be chosen to calculate the vulnerability probability of gas pipelines in karst areas. Firstly, the reliability index β is introduced. Then, it is assumed that the vulnerability function V obeys a normal distribution, with an average value of μv and a standard deviation of σv. In this regard, the vulnerability probability function is shown in Formula (17). (17) If v = μv + σvu, then dv = σvdu. When v = 0, then u = -μv/σv. When v → -∞, then u → -∞. In this regard, the vulnerability probability function can be further transformed into Formula (18). (18) Where Φ(β) is the cumulative function of standard normal distribution; β is the pipeline reliability index, and β = μv/σv. For the general form of vulnerability function: (19) Where v1, v2, …, vn are n independent normal random variables, with an average of μj (j = 1, 2, …, n) and a standard deviation of σj (j = 1, 2, …, n). Assuming that a design checkpoint is represented as v* = (v1*, v2*, …, vn*), and g(v*) is equal to 0. Then, the functional function is expanded according to the Taylor series, and its first order expansion is shown in Formula (20). (20) The mean and standard deviation of ML are shown in Formula (21) and Formula (22). (21)(22) On this basis, the expression of reliability index β is shown in Formula (23). (23) The disaster-causing possibility of gas pipelines in karst areas is comprehensively investigated from two dimensions: the hazard of disaster events and the vulnerability of gas pipelines. If the disaster-causing probability is P, then P = P(D)P(V). According to the standard "API 581–2008 Risk-based Inspection Technology" issued by the American Petroleum Institute (API) and "DNV-RP-F116 Integrity Management of Submarine Pipeline Systems" issued by the Det Norske Veritas (DNV), the disaster-causing probability of gas pipelines in karst areas is divided into five grades, as shown in Table 6. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 6. Classification criteria for disaster-causing probability. https://doi.org/10.1371/journal.pone.0316820.t006 3. Case analysis Zhijin County, Guizhou Province, China is located at the intersection of the Nayong-Kaiyang east-west structural belt and the Zhijin northeast structural belt. Under the joint sculpture of internal and external forces, various karst landscapes have been formed in the region. Under the influence of multiple factors such as disorderly mining, underground cave development, sustained heavy rainfall, and seismic activity, karst collapse occurs frequently. A buried gas pipeline project in Zhijin county is selected as the study case. 3.1 Calculation of hazard probability In this study, six experts from different fields such as gas pipeline design, construction, management, safety risk assessment, disaster rescue, disaster prevention and control are hired to form an expert evaluation team. The descriptive statistics of the experts are shown in Table 7. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 7. Descriptive statistics of the experts. https://doi.org/10.1371/journal.pone.0316820.t007 Then, the optimization matrixes of knowledge level, experience level, information source and justice degree of different experts are constructed (Data in S1 File). Through the hierarchical analysis of expert authority evaluation, the evaluation ability weights of each expert are obtained as w = (we1, we2, we3, we4, we5, we6) = (0.1720, 0.1422, 0.1715, 0.1928, 0.1441, 0.1773). Then, all the experts are hired to evaluate the hazard level for various basic events based on the knowledge and experience in their respective fields of work (Data in S2 File). On this basis, the fuzzy evaluation language of the experts is quantified based on the membership function, so as to investigate the fuzzy probability of various basic events. In addition, considering that the events related to karst collapse are in a critical position in the hazard evolution of gas pipeline disasters, the risk level of basic events t15~t19 is classified combined with the actual exploration data in karst areas. Meanwhile, micro concretization is carried out for each basic event, and the overall evaluation results are obtained through weighted average processing. Through this approach, it is helpful to more accurately evaluate the influence of key disaster events on hazard probability combined with engineering field data. The specific evaluation criteria of events related to karst collapse are shown in Table 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 8. The specific evaluation criteria of events related to karst collapse. https://doi.org/10.1371/journal.pone.0316820.t008 Taking "t1 Low level of education" as an example, the evaluation results of the expert evaluation team on the hazard level of this basic event are L, L, RL, RL, L and RL respectively. Based on the Formula (5), the membership functions of L and RL are constructed as shown in Formula (24) and Formula (25) respectively. On this basis, the corresponding α-cut sets are EWα = [0.15+0.1α, 0.45–0.1α] and ESα = [0.55+0.1α, 0.85–0.1α] respectively. (24)(25) The fuzzy number M considering the expert evaluation authority is further obtained as shown in Formula (26). (26) Combined with fuzzy set extension theory, the comprehensive evaluation function of the fuzzy number M considering the weights of different experts is shown in Formula (27). (27) Based on the Formula (6) to Formula (10), the left fuzzy possibility value FSP,L(M) is calculated to be 0.8715, and the right fuzzy possibility value FSP,R(M) is calculated to be 0.3258. Combining the left and right fuzzy possibility values, the fuzzy possibility score is calculated to be 0.2271. Based on the Formula (11) and Formula (12), the prior probability value pj(1) of the basic event "t1 Low level of education" is calculated to be 0.0003. Similarly, the prior probabilities of other disaster events are calculated. Then, according to the Bayesian network calculation rule, the probability distribution of nodes at all levels is calculated step by step. The specific results are summarized in Fig 9. Finally, the hazard probability P(D) of gas pipeline disaster events in the study area is calculated to be 0.0352. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Calculation results of hazard probability based on the Bayesian network model. https://doi.org/10.1371/journal.pone.0316820.g009 3.2 Calculation of vulnerability probability In order to determine the critical value of each vulnerability index, the finite element model of buried pipelines in karst areas is established based on the ABAQUS software, as shown in Fig 10. Considering that gas pipelines mainly adopt thin-walled structures, the 4-node shell elements are used for pipeline modeling to simulate their actual characteristics. In order to ensure the simulation accuracy and control the calculation cost, the pipeline finite element model is divided into 100 units in the longitudinal direction and 30 units in the circumferential direction, totaling 3000 shell units. The basic parameters of PE80 pipeline are shown in Table 9. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. Finite element model of buried pipeline in karst areas. https://doi.org/10.1371/journal.pone.0316820.g010 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 9. The basic parameters of PE80 pipeline. https://doi.org/10.1371/journal.pone.0316820.t009 When the soil is in the elastic-plastic stage, it will show the nonlinear characteristics of materials under the multiple influences of external environment, gravity and pipeline reaction force. Therefore, the Mohr-Coulomb model is chosen to simulate soil properties. In addition, the 8-node 3D solid elements are used for soil modeling. They have the characteristics of elastoplasticity and large displacement, and conforms to the characteristics of karst collapse. According to the study on design parameters by the American Society of Mechanical Engineers [34], the effective calculation length, width, and height of the soil model are set to 20m, 16m, and 12m, respectively. Meanwhile, the structured grid division technology is used to mesh the soil elements within the vertical and horizontal 1m of the pipeline, and finally the soil model is divided into 127,008 solid elements. The property parameters of the soil model are based on the field exploration data from karst areas in Guizhou Province, China. Considering the complex and variable geological structure of karst areas, the soil layers are simplified into three categories: plain fill, clay, and bedrock. The specific parameters are shown in Table 10. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 10. The property parameters of soil model. https://doi.org/10.1371/journal.pone.0316820.t010 In the process of karst collapse, the collapse range of the upper soil layer gradually expands with the development of roof cracks. That is, the supporting force on the karst overburden near the crack gradually decreases in the vertical direction, rather than the whole soil directly moving downwards. The element birth and death technique provided by the ABAQUS finite element analysis software can achieve the "killing" and reactivation of soil elements, thereby accurately simulating the soil loss during the gradual development of karst collapse. Therefore, the karst soil in the lower part of the pipeline is divided into five unit sets from inside to outside along the collapse direction, and then the "Invalid in this step" is set in the software interaction attribute module. On this basis, the gradual development characteristics of karst collapse from the central soil cave to both sides are simulated based on the element birth and death technique. The mechanical characteristics of buried gas pipeline under the influence of karst collapse are shown in Fig 11. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. The mechanical characteristics of buried gas pipeline under the influence of karst collapse. https://doi.org/10.1371/journal.pone.0316820.g011 Based on the finite element simulation results, the value ranges of different vulnerability parameters are further analyzed: (1) Karst collapse parameters When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, a wall thickness of 18.2mm and a buried depth of 1m, the critical length of karst soil collapse that caused the pipeline failure is 5.6m, the critical width of karst soil collapse is 6.1m, and the critical thickness of karst overburden is 2.6m. At this time, for Formula (16), a1 = Kscl/5.6, a2 = Kscw/6.1 and a3 = Kot/2.6. Combined with the karst exploration data along the gas pipelines in the study area, the length of karst soil cave ranges from 2 to 10m, the width of karst soil cave ranges from 0.3 to 6.5m, and the thickness of karst overburden ranges from 1.8 to 19.3m. In this regard, the value ranges of a1, a2 and a3 can be obtained. It is assumed that the a1, a2 and a3 are randomly distributed in their respective value ranges, and 10000 random numbers are generated in their value ranges respectively based on the Matlab programming (Data in S3 File). The average values and standard deviations of each group of random numbers are calculated and summarized in Table 11. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 11. Related values of vulnerability parameters. https://doi.org/10.1371/journal.pone.0316820.t011 (2) Pipeline operation parameters When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, and a buried depth of 1m, the critical pipe wall thickness is 16.89mm. At this time, for Formula (16), a4 = Pwt/16.88. According to the national standard of the People’s Republic of China, "GB 15558.1 Buried Polyethylene (PE) Pipeline Systems for Gas Use Part 1: Pipes", the minimum wall thickness of SDR11 series pipeline materials is 18.2mm, and the deviation between the wall thickness at any point and the minimum wall thickness is within 2mm. Therefore, the pipe wall thickness ranges from 16.2 to 20.2mm. In this regard, the value range of the a4 can be obtained. It is assumed that the a4 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm and a wall thickness of 18.2mm, the critical pipe buried depth that caused the pipeline failure is 0.97m. At this time, for Formula (16), a5 = Pbd/0.97. Due to the complex and changeable geological environment along the gas pipeline in the study area, the buried depth of some pipeline sections exceeds or is less than the standard depth by about 1m during the process of actual laying. On the whole, the buried depth is in the range of 0.6 ~ 1.5m. In this regard, the value range of the a5 can be obtained. It is assumed that the a5 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. According to the national standard of the People’s Republic of China, "GB 50494–2009 Technical Specification for Urban Gas", the design service life of gas pipelines in China should be greater than 30 years, and the service life of PE pipelines can generally reach 40 ~ 50 years. Therefore, the service life of PE gas pipelines is taken as 50 years in this study. At this time, for Formula (16), a8 = Lps/50. The PE gas pipeline in the study area has been in service for 4 years. In this regard, the value range of the a8 can be obtained. It is assumed that the a8 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. (3) Soil property parameters When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, a wall thickness of 18.2mm and a buried depth of 1m, the critical soil elastic modulus is 22.63MPa. At this time, for Formula (16), a6 = Sem/22.63. Due to the complex geological structure along the gas pipeline in the study area, the soil layer where the buried pipeline is located is mainly composed of plain fill and clay, and its elastic modulus is mainly in the range of 15 ~ 25MPa. In this regard, the value range of the a6 can be obtained. It is assumed that the a6 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, a wall thickness of 18.2mm and a buried depth of 1m, the critical friction coefficient between pipeline and soil is 0.45. At this time, for Formula (16), a7 = Fc/0.45. Due to the complex and variable soil properties and moisture content along the gas pipeline in the study area, the friction coefficient is mainly in the range of 0.3 ~ 0.5. In this regard, the value range of the a7 can be obtained. It is assumed that the a7 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. The initial design checkpoint is set based on the average value of each vulnerability parameters, that is, v*(0) = (1.07, 0.55, 0.79, 0.94, 1.13, 0.54, 0.88, 0.89). Based on the Formula (17) to Formula (23), the corresponding calculation program is compiled based on the MATLAB programming. The iteration number of model calculation is 13, and the reliability index β is finally calculated to be -1.2639. By consulting the standard normal distribution table, the vulnerability probability P(V) of gas pipelines in karst areas based on the design checkpoint method is 0.1038, as shown in Fig 12 [35]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Relationship between the vulnerability probability P(V) and the reliability index β. https://doi.org/10.1371/journal.pone.0316820.g012 In this study, the disaster-causing possibility of gas pipelines in karst areas is comprehensively measured from two dimensions: the hazard of disaster events and the vulnerability of gas pipelines. Therefore, the disaster-causing probability is as follows: P = P(D)P(V) = 0.0352 × 0.1038 = 0.0037. According to the classification criteria, the disaster-causing probability is at level III (It happens occasionally). This level is unacceptable for the daily operation of gas pipelines, so it is necessary to take corresponding disaster prevention and control measures to reduce the risk of disaster events and the vulnerability of gas pipelines in karst areas. 3.1 Calculation of hazard probability In this study, six experts from different fields such as gas pipeline design, construction, management, safety risk assessment, disaster rescue, disaster prevention and control are hired to form an expert evaluation team. The descriptive statistics of the experts are shown in Table 7. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 7. Descriptive statistics of the experts. https://doi.org/10.1371/journal.pone.0316820.t007 Then, the optimization matrixes of knowledge level, experience level, information source and justice degree of different experts are constructed (Data in S1 File). Through the hierarchical analysis of expert authority evaluation, the evaluation ability weights of each expert are obtained as w = (we1, we2, we3, we4, we5, we6) = (0.1720, 0.1422, 0.1715, 0.1928, 0.1441, 0.1773). Then, all the experts are hired to evaluate the hazard level for various basic events based on the knowledge and experience in their respective fields of work (Data in S2 File). On this basis, the fuzzy evaluation language of the experts is quantified based on the membership function, so as to investigate the fuzzy probability of various basic events. In addition, considering that the events related to karst collapse are in a critical position in the hazard evolution of gas pipeline disasters, the risk level of basic events t15~t19 is classified combined with the actual exploration data in karst areas. Meanwhile, micro concretization is carried out for each basic event, and the overall evaluation results are obtained through weighted average processing. Through this approach, it is helpful to more accurately evaluate the influence of key disaster events on hazard probability combined with engineering field data. The specific evaluation criteria of events related to karst collapse are shown in Table 8. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 8. The specific evaluation criteria of events related to karst collapse. https://doi.org/10.1371/journal.pone.0316820.t008 Taking "t1 Low level of education" as an example, the evaluation results of the expert evaluation team on the hazard level of this basic event are L, L, RL, RL, L and RL respectively. Based on the Formula (5), the membership functions of L and RL are constructed as shown in Formula (24) and Formula (25) respectively. On this basis, the corresponding α-cut sets are EWα = [0.15+0.1α, 0.45–0.1α] and ESα = [0.55+0.1α, 0.85–0.1α] respectively. (24)(25) The fuzzy number M considering the expert evaluation authority is further obtained as shown in Formula (26). (26) Combined with fuzzy set extension theory, the comprehensive evaluation function of the fuzzy number M considering the weights of different experts is shown in Formula (27). (27) Based on the Formula (6) to Formula (10), the left fuzzy possibility value FSP,L(M) is calculated to be 0.8715, and the right fuzzy possibility value FSP,R(M) is calculated to be 0.3258. Combining the left and right fuzzy possibility values, the fuzzy possibility score is calculated to be 0.2271. Based on the Formula (11) and Formula (12), the prior probability value pj(1) of the basic event "t1 Low level of education" is calculated to be 0.0003. Similarly, the prior probabilities of other disaster events are calculated. Then, according to the Bayesian network calculation rule, the probability distribution of nodes at all levels is calculated step by step. The specific results are summarized in Fig 9. Finally, the hazard probability P(D) of gas pipeline disaster events in the study area is calculated to be 0.0352. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 9. Calculation results of hazard probability based on the Bayesian network model. https://doi.org/10.1371/journal.pone.0316820.g009 3.2 Calculation of vulnerability probability In order to determine the critical value of each vulnerability index, the finite element model of buried pipelines in karst areas is established based on the ABAQUS software, as shown in Fig 10. Considering that gas pipelines mainly adopt thin-walled structures, the 4-node shell elements are used for pipeline modeling to simulate their actual characteristics. In order to ensure the simulation accuracy and control the calculation cost, the pipeline finite element model is divided into 100 units in the longitudinal direction and 30 units in the circumferential direction, totaling 3000 shell units. The basic parameters of PE80 pipeline are shown in Table 9. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 10. Finite element model of buried pipeline in karst areas. https://doi.org/10.1371/journal.pone.0316820.g010 Download: PPT PowerPoint slide PNG larger image TIFF original image Table 9. The basic parameters of PE80 pipeline. https://doi.org/10.1371/journal.pone.0316820.t009 When the soil is in the elastic-plastic stage, it will show the nonlinear characteristics of materials under the multiple influences of external environment, gravity and pipeline reaction force. Therefore, the Mohr-Coulomb model is chosen to simulate soil properties. In addition, the 8-node 3D solid elements are used for soil modeling. They have the characteristics of elastoplasticity and large displacement, and conforms to the characteristics of karst collapse. According to the study on design parameters by the American Society of Mechanical Engineers [34], the effective calculation length, width, and height of the soil model are set to 20m, 16m, and 12m, respectively. Meanwhile, the structured grid division technology is used to mesh the soil elements within the vertical and horizontal 1m of the pipeline, and finally the soil model is divided into 127,008 solid elements. The property parameters of the soil model are based on the field exploration data from karst areas in Guizhou Province, China. Considering the complex and variable geological structure of karst areas, the soil layers are simplified into three categories: plain fill, clay, and bedrock. The specific parameters are shown in Table 10. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 10. The property parameters of soil model. https://doi.org/10.1371/journal.pone.0316820.t010 In the process of karst collapse, the collapse range of the upper soil layer gradually expands with the development of roof cracks. That is, the supporting force on the karst overburden near the crack gradually decreases in the vertical direction, rather than the whole soil directly moving downwards. The element birth and death technique provided by the ABAQUS finite element analysis software can achieve the "killing" and reactivation of soil elements, thereby accurately simulating the soil loss during the gradual development of karst collapse. Therefore, the karst soil in the lower part of the pipeline is divided into five unit sets from inside to outside along the collapse direction, and then the "Invalid in this step" is set in the software interaction attribute module. On this basis, the gradual development characteristics of karst collapse from the central soil cave to both sides are simulated based on the element birth and death technique. The mechanical characteristics of buried gas pipeline under the influence of karst collapse are shown in Fig 11. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 11. The mechanical characteristics of buried gas pipeline under the influence of karst collapse. https://doi.org/10.1371/journal.pone.0316820.g011 Based on the finite element simulation results, the value ranges of different vulnerability parameters are further analyzed: (1) Karst collapse parameters When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, a wall thickness of 18.2mm and a buried depth of 1m, the critical length of karst soil collapse that caused the pipeline failure is 5.6m, the critical width of karst soil collapse is 6.1m, and the critical thickness of karst overburden is 2.6m. At this time, for Formula (16), a1 = Kscl/5.6, a2 = Kscw/6.1 and a3 = Kot/2.6. Combined with the karst exploration data along the gas pipelines in the study area, the length of karst soil cave ranges from 2 to 10m, the width of karst soil cave ranges from 0.3 to 6.5m, and the thickness of karst overburden ranges from 1.8 to 19.3m. In this regard, the value ranges of a1, a2 and a3 can be obtained. It is assumed that the a1, a2 and a3 are randomly distributed in their respective value ranges, and 10000 random numbers are generated in their value ranges respectively based on the Matlab programming (Data in S3 File). The average values and standard deviations of each group of random numbers are calculated and summarized in Table 11. Download: PPT PowerPoint slide PNG larger image TIFF original image Table 11. Related values of vulnerability parameters. https://doi.org/10.1371/journal.pone.0316820.t011 (2) Pipeline operation parameters When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, and a buried depth of 1m, the critical pipe wall thickness is 16.89mm. At this time, for Formula (16), a4 = Pwt/16.88. According to the national standard of the People’s Republic of China, "GB 15558.1 Buried Polyethylene (PE) Pipeline Systems for Gas Use Part 1: Pipes", the minimum wall thickness of SDR11 series pipeline materials is 18.2mm, and the deviation between the wall thickness at any point and the minimum wall thickness is within 2mm. Therefore, the pipe wall thickness ranges from 16.2 to 20.2mm. In this regard, the value range of the a4 can be obtained. It is assumed that the a4 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm and a wall thickness of 18.2mm, the critical pipe buried depth that caused the pipeline failure is 0.97m. At this time, for Formula (16), a5 = Pbd/0.97. Due to the complex and changeable geological environment along the gas pipeline in the study area, the buried depth of some pipeline sections exceeds or is less than the standard depth by about 1m during the process of actual laying. On the whole, the buried depth is in the range of 0.6 ~ 1.5m. In this regard, the value range of the a5 can be obtained. It is assumed that the a5 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. According to the national standard of the People’s Republic of China, "GB 50494–2009 Technical Specification for Urban Gas", the design service life of gas pipelines in China should be greater than 30 years, and the service life of PE pipelines can generally reach 40 ~ 50 years. Therefore, the service life of PE gas pipelines is taken as 50 years in this study. At this time, for Formula (16), a8 = Lps/50. The PE gas pipeline in the study area has been in service for 4 years. In this regard, the value range of the a8 can be obtained. It is assumed that the a8 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. (3) Soil property parameters When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, a wall thickness of 18.2mm and a buried depth of 1m, the critical soil elastic modulus is 22.63MPa. At this time, for Formula (16), a6 = Sem/22.63. Due to the complex geological structure along the gas pipeline in the study area, the soil layer where the buried pipeline is located is mainly composed of plain fill and clay, and its elastic modulus is mainly in the range of 15 ~ 25MPa. In this regard, the value range of the a6 can be obtained. It is assumed that the a6 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. When the karst collapse load directly affects the PE80 pipeline with a diameter of 200mm, a wall thickness of 18.2mm and a buried depth of 1m, the critical friction coefficient between pipeline and soil is 0.45. At this time, for Formula (16), a7 = Fc/0.45. Due to the complex and variable soil properties and moisture content along the gas pipeline in the study area, the friction coefficient is mainly in the range of 0.3 ~ 0.5. In this regard, the value range of the a7 can be obtained. It is assumed that the a7 is randomly distributed in its respective value range, and 10000 random numbers are generated in its value range based on the Matlab programming (Data in S3 File). The average value and standard deviation of the group of random numbers are calculated and summarized in Table 11. The initial design checkpoint is set based on the average value of each vulnerability parameters, that is, v*(0) = (1.07, 0.55, 0.79, 0.94, 1.13, 0.54, 0.88, 0.89). Based on the Formula (17) to Formula (23), the corresponding calculation program is compiled based on the MATLAB programming. The iteration number of model calculation is 13, and the reliability index β is finally calculated to be -1.2639. By consulting the standard normal distribution table, the vulnerability probability P(V) of gas pipelines in karst areas based on the design checkpoint method is 0.1038, as shown in Fig 12 [35]. Download: PPT PowerPoint slide PNG larger image TIFF original image Fig 12. Relationship between the vulnerability probability P(V) and the reliability index β. https://doi.org/10.1371/journal.pone.0316820.g012 In this study, the disaster-causing possibility of gas pipelines in karst areas is comprehensively measured from two dimensions: the hazard of disaster events and the vulnerability of gas pipelines. Therefore, the disaster-causing probability is as follows: P = P(D)P(V) = 0.0352 × 0.1038 = 0.0037. According to the classification criteria, the disaster-causing probability is at level III (It happens occasionally). This level is unacceptable for the daily operation of gas pipelines, so it is necessary to take corresponding disaster prevention and control measures to reduce the risk of disaster events and the vulnerability of gas pipelines in karst areas. 4. Conclusion (1) Based on the disaster system theory, the hazard evaluation index system of gas pipeline disaster events in karst areas is established from three aspects: the activity of disaster-prone environment, the risk of disaster factors and the vulnerability of disaster-bearing bodies. Meanwhile, based on the Bayesian network model, the calculation of the occurrence probability of basic events is equivalent to the propagation and renewal of beliefs, which can make up for the limitations of effective calculation bits and complicated calculation process of fault tree model. Therefore, it is helpful to measure the hazard probability of disaster events more simply and accurately from a systematic perspective. (2) Considering the pipeline vulnerability level and its relationship with the disaster type and intensity, pipeline structure and function, disaster and pipeline space-time configuration, a mathematical model system is established based on the polyethylene pipeline characteristics. Combined with the design checkpoint method, the vulnerability probability of gas pipelines is calculated from the perspective of the interaction between disaster hazard intensity and pipeline resistance capacity. It can make up for the subjective deficiency in previous research, which mainly relied on expert experience for evaluation. (3) Based on the data of karst exploration and field investigation, combined with fuzzy theory, the probability distribution of nodes at different levels of the Bayesian network is calculated step by step, and the hazard probability of disaster events is calculated to be 0.0352. Based on the finite element simulation results, combined with the reliability analysis theory, the vulnerability probability of gas pipelines is calculated to be 0.1038. The disaster-causing possibility of gas pipelines in karst areas is comprehensively measured from two dimensions: the hazard of disaster events and the vulnerability of gas pipelines, and finally it is calculated to be 0.0037. Supporting information S1 File. The optimization matrixes of knowledge level, experience level, information source and justice degree. https://doi.org/10.1371/journal.pone.0316820.s001 (PDF) S2 File. Expert evaluation results of hazard level. https://doi.org/10.1371/journal.pone.0316820.s002 (PDF) S3 File. Random numbers of 8 vulnerable parameters. https://doi.org/10.1371/journal.pone.0316820.s003 (XLSX) TI - Study on disaster-causing probability evaluation of gas pipeline in karst area JF - PLoS ONE DO - 10.1371/journal.pone.0316820 DA - 2025-02-03 UR - https://www.deepdyve.com/lp/public-library-of-science-plos-journal/study-on-disaster-causing-probability-evaluation-of-gas-pipeline-in-6SKgpgAMAT SP - e0316820 VL - 20 IS - 2 DP - DeepDyve ER -