Real-time Risk Assessment for Petrochemical Industry

Prévia do material em texto

Reliability Engineering and System Safety 226 (2022) 108700
Available online 2 July 2022
0951-8320/© 2022 Elsevier Ltd. All rights reserved.
A real-time probabilistic risk assessment method for the petrochemical 
industry based on data monitoring 
Rui He a,*, Jingyu Zhu b, Guoming Chen b, Zhigang Tian a 
a Department of Mechanical Engineering, University of Alberta, Edmonton, AB, Canada 
b Centre for Offshore Engineering and Safety Technology (COEST), China University of Petroleum (East China), No.66, Changjiang West Road, Qingdao, China 
A R T I C L E I N F O 
Keywords: 
Probabilistic risk assessment (PRA) 
Dynamic Bayesian network (DBN) 
Data monitoring 
Petrochemical industry 
A B S T R A C T 
Safety is of high societal concern in the petrochemical industry. With advancing digitalization, the techniques of 
probabilistic risk assessment (PRA), which provide a system-level perspective for industrial risk analysis, have 
become increasingly dynamic. This paper provides a further advancement in this direction by proposing a risk 
updating method based on the dynamic Bayesian network (DBN) to incorporate data monitoring into PRA in real- 
time. We update the probabilistic risk for basic events in Bayesian networks based on their prior probability and 
the probability of their online monitoring signals exceeding the alarm threshold. The key idea behind this 
approach is that if the signal of a basic event exceeds the alarm threshold, the corresponding risk should increase. 
Otherwise, the basic event should have a low risk. The contribution in this paper is twofold. First, residual 
methods are developed to estimate signal probabilities exceeding abnormal thresholds. Second, a DBN model is 
proposed to integrate prior risk with data monitoring for risk analysis. The proposed DBN model does not require 
additional expert knowledge and historical accident data to define the conditional relationship between data 
monitoring and prior risk. The proposed approach is validated using the RT 580 experimental setup and managed 
pressure drilling operations. 
1. Introduction 
Modern petrochemical processes are complex socio-technical sys-
tems that are susceptible to catastrophic accidents due to an abnormal 
process, instrumentation, or human error [1]. Probabilistic risk assess-
ment (PRA) provides a well-established family of techniques to address 
safety, risk, and reliability of these complex petrochemical systems [2]. 
To equip PRA with the capability of investigating accident propagation 
mechanisms, various methodologies have been proposed over the past 
decades, including fault tree analysis (FTA) [3], event tree analysis 
(ETA) [4], bow-tie (BT) [5], and Bayesian network (BN) [6]. These 
traditional PRA techniques consider possible accident propagation with 
chain rules and use pre-set deterministic failure probabilities to assess 
the likelihood of accidents. However, traditional PRA methods are 
limited by the use of static failure probabilities and thus is hard to 
capture deviations and changes in system status. As the modern industry 
moves towards a data-driven “Industry 4.0”, advances in sensor tech-
nology have made PRA techniques more and more dynamic, but effec-
tively utilizing data information in PRA remains a challenge [7]. 
Most petrochemical processes are intermittent systems with multiple 
working phases. During system operations, the risk of accidents may 
change over time, and the accident propagation mechanism during 
phase shifts may be different. To account for the time-dependent 
behavior of accidents, engineers started to consider additional time 
slices in PRA methods. OREDA [8] provides a reference for exponential 
failure rates of most petrochemical equipment, and thus it is possible to 
develop static PRA into dynamic forms based on failure rates. Typically, 
constant failure rates have been widely used for extending BN into a 
dynamic Bayesian network (DBN) to consider time-dependent opera-
tions and equipment degradations [9–11]. When statistical accident 
data or precursor data is available, the risk can also be updated in PRA 
by Bayesian inference. Exponential distribution and Beta distribution 
can be assumed as the basic probability distribution of possible accidents 
in ETA or BN, and the failure probabilities of safety barriers are updated 
by counting the number of accidents [12,13]. Since the available acci-
dent data is often limited, hierarchical Bayesian analysis (HBA) is pro-
posed to use scarce statistical data to update the failure probability 
distribution of safety barriers [14,15]. Other graph-based methods, such 
as directed acyclic graphs (DAG) that assign a deterministic value to 
each path for the risk propagation analysis also use statistical accident 
* Corresponding author. 
E-mail address: rhe5@ualberta.ca (R. He). 
Contents lists available at ScienceDirect 
Reliability Engineering and System Safety 
journal homepage: www.elsevier.com/locate/ress 
https://doi.org/10.1016/j.ress.2022.108700 
Received 20 August 2021; Received in revised form 18 May 2022; Accepted 26 June 2022 
mailto:rhe5@ualberta.ca
www.sciencedirect.com/science/journal/09518320
https://www.elsevier.com/locate/ress
https://doi.org/10.1016/j.ress.2022.108700
https://doi.org/10.1016/j.ress.2022.108700
https://doi.org/10.1016/j.ress.2022.108700
http://crossmark.crossref.org/dialog/?doi=10.1016/j.ress.2022.108700&domain=pdf
Reliability Engineering and System Safety 226 (2022) 108700
2
data to update risk by assuming distributions of risk factors [16]. One 
issue is that statistical data are often collected from a population of 
similar systems and cannot fully interpret the specific features of target 
systems. In addition, the accident data is commonly hard to collect in the 
short term, and thus these PRA techniques can only be updated regularly 
over long periods, such as yearly or weekly [17]. 
Due to the real-time characteristics of operating parameters in 
petrochemical industries, the credibility of statistical data is always in 
question [18]. Advances in sensor technology provide online data sup-
port for PRA techniques, and therefore, PRA methods have begun to 
update estimated risk in real-time using online measurements. Unlike 
process monitoring, PRA uses risk to measure the lack of safety from 
system-level perspectives and thus can identify possible problems and 
make corrective decisions before accidents occur [19]. Commonly, an 
accident is caused by a combination of several basic events, including 
abnormal operating conditions and safety barrier failures. Condition 
data can be incorporated in degradation models based on reliability 
assessment methods when considering the dynamic mechanism of safety 
barriers. Initially, a condition fault tree is proposed to update the failure 
rates of components according to equipment condition data [20]. 
Model-based methods, like Kalman filter and particle filtering (PF), are 
developed to update the degradation state of safety barriers based on 
real-time measurements [21–23]. In addition, Gaussian mixture model 
(GMM) is developed in [24] to estimate the degradation of safety bar-
riers and identify health status in real-time. Furthermore, in [7], a sys-
tematical framework is presented to incorporate the degradation model 
into PRA techniques. Even though the degradation characteristics of 
safety barriers are effectively modeled, the reliability-based methods 
cannot fully meet the requirement of petrochemical industries because 
the instantaneous changes in operations cannot be considered. In addi-
tion, the process equipment will not be significantly degraded in a short 
operating period. 
To equip PRA with data monitoring, some works construct data- 
driven risk models based on large amounts of historical or experi-
mental accident data. Commonly, as in [25–27], BN structural and 
parameter learning are performed to learn probabilistic relationships 
between operating parameters andaccidents using historical event data, 
and the learned BN can be used to predict the probability of accidents 
when the online data is available. In [28–30], data monitoring is 
incorporated into BN models for risk analysis, where BN parameters are 
determined using statistical methods such as maximum likelihood esti-
mation (MLE) based on sufficient accident-related data. Also in [31,32], 
a risk-based warning model is constructed for real-time risk analysis 
using historical data and expert knowledge. Oftentimes, 
accident-related data are challenging to collect, and thus statistical and 
data-driven approaches may lead to biased and inconsistent estimates in 
PRA models. Other related works determine different failure probabili-
ties for basic events under varying operations. In [33,34], different 
risk-based loss functions are defined to use monitoring data to assess 
accident probabilities in real-time. In addition, in [35–37], the online 
monitoring signal is divided into multiple discrete levels (e.g., high, 
medium, low), and corresponding failure probabilities are set for basic 
events under different groups according to expert knowledge. However, 
despite no use of historical data, these PRA techniques determine the 
risks of basic events based only on the online information while ignoring 
their prior probabilities in real-time risk estimations. If the prior prob-
ability is not considered, the risk of an event will depend entirely on the 
data monitoring and will be unrelated to its own accident mechanism. 
Only a few works, for example [38,39], consider both data monitoring 
and prior probabilities in PRA without using historical accident data, but 
expert knowledge is still required to determine the effect of different 
monitoring parameters on the probability of basic events or accidents. 
In this work, we provide a further advancement in this direction by 
proposing a real-time PRA approach that embeds both data monitoring 
and prior knowledge in risk analysis. We start with the estimates of the 
probability that online signals exceed the alarm threshold. Thereafter, 
these estimated probabilities are incorporated into PRA to update the 
risk of basic events based on their prior knowledge. The significant 
contribution of this work is twofold. First, residual-based methods are 
developed for environmental and process monitoring to estimate the 
probability of their online signals exceeding abnormal thresholds. Sec-
ond, a DBN model is proposed to integrate prior risk with data moni-
toring for risk analysis. The proposed DBN model does not require 
additional expert knowledge and historical accident data for each basic 
event to determine the conditional relationship (conditional probability 
table) between data monitoring and prior risk so that the proposed PRA 
is practical, especially when historical data is scarce. 
The rest of the paper is organized as follows. In Section 2, the pre-
liminaries and engineering motivations are briefly introduced, and the 
critical problems are formulated. In Section 3, the proposed real-time 
PRA method is described. Section 4 uses two cases, including the RT 
580 experimental setup and managed pressure drilling (MPD) operation, 
to demonstrate and validate the proposed method. Finally, conclusions 
are made in Section 5. 
2. Preliminaries and problem formulation 
This work aims to predict the probability of accidents and to search 
for critical events with potential accident risks. Let TE represents the top 
event (main accident) of the Fault Tree (FT) and suppose that there are N 
basic events to determine the accident, denoted by BEi, i = 1, 2, …, N. 
Based on the logical relationship of FT gates, the probability of an ac-
cident is a function of the prior probability of basic events 
PTE = fFT(PBE1 ,PBE2 ,…,PBEN ), (1) 
wherefFT(⋅) denotes the FT model and PBEi is the prior probability of the 
i-th basic event obtained by FTA. 
The mapping algorithm in [6] can be used to set conditional prob-
ability tables (CPTs) to map FT to BN. In order to discover the potential 
basic events that may lead to accidents, the risk index is defined as the 
conditional probability of BEi given that TE has already occurred. The 
conditional probabilities can be estimated based on the rules of BN 
P(BEi|TE) = P{BEioccurs|TEhasoccured}
=
P(TE|BEi)PBEi∑
BEi
P(TE,BEi)
, i = 1, 2,…N. (2) 
To meet the timing requirement of petrochemical industries, data 
monitoring is considered in PRA. Fig. 1 shows a simple process (feeding 
tank) and its FT with high-level incident [40]. The top event, i.e., high 
level, results from the combined effect of both the abnormal operating 
flow and the failure of safety barriers. In this case, the condition of input 
and output flow and the level control logic can be respectively observed 
from FI01, FI02, and LIC01, and the pump status can be assessed based 
on the accelerometers. The online-monitoring data can reflect the sys-
tem operating status. For example, if FI01 is high, the risk of increased 
feed should also rise. However, since different events have different 
prior probabilities, the updated risk should not only depend on data 
monitoring but also on the prior risk of the event. 
When a BN model is constructed, the nodes of BN are categorized 
into static nodes and dynamic nodes. The real-time probability of a static 
node is always equal to its prior probability, while the probability of 
dynamic nodes will be continuously updated based on data monitoring. 
Dynamic nodes can be further classified into environmental nodes and 
process nodes according to signal characteristics. Environmental nodes 
refer to various environmental variables that are not affected by changes 
in operating conditions, such as wind speed and formation pressure. 
Process nodes, such as input flow, tank level, and pump vibration, are 
associated with process variables affected by varying working condi-
tions. We assume that the risk of m environmental nodes can be updated 
based on environmental monitoring, and the risk of n process nodes can 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
3
be updated based on process monitoring. In this work, the probability of 
online signals beyond the alarm threshold is first estimated for each 
dynamic node, which can be regarded as the probability that the event 
deviates from the normal status, denoted byPMONi , i = 1, 2,…, m+n. 
Afterward, the risks of dynamic nodes are estimated in real-time by 
integrating PMON and their prior probability PBE. The PRA tasks are 
formally defined as: 
(1) forward risk updating: update the probability of basic events at 
the current time t = tk, k = 1, 2, …,T, based on data monitoring 
available up to tk, to predict the probability of top event at time tk; 
(2) posterior risk estimation: given that the top event has occurred, 
estimate the risk index (posterior probability) of basic events at 
time tk to discover the critical event with potential accident risks. 
Fig. 1. Process diagram and Fault Tree model for a feeding tank. 
Fig. 2. Procedures for real-time PRA based on data monitoring 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
4
3. Proposed methodology 
Fig. 2 summarizes the major steps of the proposed real-time PRA 
method for petrochemical industries. Before conducting a risk analysis, 
the event type of each basic event is first determined. After that, the 
probability of the online signal exceeding the alarm thresholdPMON is 
estimated for each dynamic node. We show how to evaluatePMONi , i = 1, 
2,…, m, for environmental nodes in Section 3.1.1, and present how to 
estimate PMONi , i = m+1, m+2,…, m+n, for process nodes in Section 
3.1.2. Finally, in Section 3.2, a novel PRA method is presented to inte-
grate PMON with prior probability PBE based on the idea of DBN for risk 
updating and posterior estimations in real-time. At time t = tk, the for-
ward risk updating is performed first, then the posterior probability of 
basic events can be estimated based on the updated probability. 
3.1. Probability calculation for online signal exceeding alarm threshold 
3.1.1. Probability estimations for environmental nodes 
As mentioned before, the environmental nodes in this section denote 
the dynamic nodes that cannot be affected by changes in system oper-
ating conditions, such as increased wind speed and deduced ambient 
temperature. This section aims to estimate the probability that the on-
line signal of the environmental node exceeds the alarm threshold. In 
general, the historical data of environment-related events is scarce for a 
new project, but the normal reference value and the alarm threshold of 
monitored parameters are often known in advance. Let AEi represent the 
i-th environmental node, i = 1, 2,…, m, and assume that there are M 
monitoring signals to determine the status of AEi denoted by xi,j, j = 1, 2, 
…, M. The probability estimation problems for an environmental node 
can be formulated as: given the reference signal value and the alarm 
threshold of xi,j, and the process operating measurements xi(t) = [xi,1(t), 
…, xi,M(t)]T, to estimate the probability of its monitoring parameters 
exceeding the alarm threshold PMONi (t)at time t. 
In order to combine multiple signals, the signal residual of the i-th 
environmental node is calculated by a multivariate logic solver 
ri,j(t) =
⃒
⃒xi,j(t) − x0,i,j
⃒
⃒
xa,i,j
, ri(t) =
1
M
∑M
j=1
ri,j(t), (3) 
where x0,i,j is the reference signal value of xi,j, xa,i,j is a positive value 
equal to the modulus of the difference between the alarm threshold and 
reference value, and therefore, ri(t) > 1 means that the residual value of 
AEi is over a threshold value, that is, xi(t) < xi,0-xi,a or xi(t) > xi,0 + xi,a. 
Let xU,i,j and xL,i,j represent the upper limit and lower limit of xi,j, 
respectively. Oftentimes, only one threshold (upper or lower limit) exists 
for a basic event. Therefore, when only xU,i,j exists and xi,j(t) > x0,i,j, or 
only xL,i,j exists and xi,j(t) < x0,i,j, the value of ri,j(t) is calculated by Eq. 
(3); otherwise, ri,j(t) is equal to 0. 
The probability PMON increases when the monitored residual moves 
further away from the normal zone. This behavior can be captured using 
a cumulative normal distribution [41]. Therefore, the probabilityPMONof 
AEi at time t is defined as 
PMONi (t) = φ
(
ri(t) − 1
σ
)
=
∫ ri(t)
− ∞
1̅̅̅
̅̅
2π
√
σ
e
(r− 1)2
2σ2 dr, (4) 
where σ = 1/3 up to the three-sigma rule. 
Fig. 3 shows a visual depiction of the probability estimation. The 
probability PMON for a point on the line ri = 0 is 0, and for a given point 
on the threshold ri = 1 is equal to 0.5. When the online residual ri(t) is at 
the threshold (ri =1), the probability of monitoring data exceeding the 
alarm threshold is 0.5 as it can either go back to the normal zone or keep 
growing and lead to the abnormal event with the same probability. 
In this way, it is possible to estimate the probability that the moni-
toring parameters of an environmental event exceed the alarm threshold 
based on data monitoring in real-time. However, PMON is not equal to the 
probability of the occurrence of environmental events. In Section 3.2, we 
combine PMON and prior knowledge to estimate the real-time risk for all 
dynamic events. 
3.1.2. Probability estimations for process nodes 
Unlike in environmental nodes, the monitoring parameters in pro-
cess nodes will be affected by varying working conditions. As operating 
conditions vary over time, the reference signal values and alarm 
thresholds of a process node may also change accordingly. Often, many 
fault-free monitoring data can be collected from systems under normal 
operating conditions. Therefore, let yi,j with j = 1, 2, …, M denotes the j- 
th process monitoring signal corresponding to the i-th process node. The 
probability estimation problems for process nodes can be formulated as: 
given historical fault-free data and online signalszi,t = [zi,1,t ,zi,2,t ,...,zi,M,t], 
i = m+1, m+2,…, m+n, to construct a nonlinear residual model for the i- 
th process node to estimate the probability of its residual signal 
exceeding the alarm threshold PMONi (t) at time t. 
Gaussian mixture model (GMM) can represent different operation 
scenarios in model building [42]. Therefore, GMM is applied to 
construct the residual model in this work. The probability density of a 
given set of training sample yi can then be approximated by GMM as 
P(yi|Λi) =
∑K
k=1
βi,kg
(
yi|λi,k
)
,
∑K
k=1
βi,k = 1, (5) 
where yi = [yi,1, …, yi,M]T is the given M-dimensional measurement 
sample, K is the number of Gaussian components, λi,k is a parameter 
vector containing mean value vector μi,k and covariance matrix Σi,k, and 
βi,k denotes the prior probability of the k-th Gaussian component. 
Conceptually, we denote the GMM compactly as Λi = {βi, μi, Σi}, where 
βi = [βi,1, …, βi,K], μi = [μi,1, …, μi,K], and Σi = [Σi,1, …, Σi,K]. 
In the offline stage, GMM parameters Λi = {βi, μi, Σi} are learned 
based on the Expectation Maximization (EM) algorithm by maximizing 
the log-likelihood of observing historical data. The minimum message 
length (MML) criterion is used to search the optimal number of Gaussian 
components. More details to train GMM can be found in [43]. 
When new signals zi,t are measured at time t in the online stage, the 
GMM model can be transformed into a dynamic form, and the proba-
bility distribution p(yi,t+1| zi,t) at the time step (t+1) can be predicted via 
the Chapman-Kolmogorov equation 
P
(
yi,t+1
⃒
⃒zi,t
)
=
∫
P
(
yi,t+1
⃒
⃒yi,t
)
P
(
yi,t
⃒
⃒zi,t
)
dyi,t
≈
∑K
k=1
βi,t,k
∫
g
(
yi,t|μi,t,k ,Σi,t,k
)
P
(
yi,t+1
⃒
⃒yi,t
)
dyi,t,
(6) 
where μtn, Σtn, and βtn denotes the updated mean value vector, 
Fig. 3. Changes in the probability of signals exceeding alarm threshold with 
monitored residual. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
5
covariance matrix, and prior probability of the n-th Gaussian component 
at the time step t, respectively. 
According to the idea of particle filter (PF) [44], the posterior 
probability distribution P(yi,t| zi,t) can be represented by a set of random 
particles y(p)i,t , p = 1, 2, …, L, with corresponding weight ω
(p)
i,t , and thus 
Eq. (6) can be represented as 
P
(
yi,t+1
⃒
⃒zi,t
)
≈
∑L
p=1
ω(p)i,t δ
(
yi,t − y(p)i,t
)
P
(
yi,t+1
⃒
⃒yi,t
)
=
∑L
p=1
ω(p)i,t P
(
yi,t+1
⃒
⃒
⃒y(p)i,t
)
,
(7) 
whereδ(⋅)stands for the Dirac delta function. 
Initially, the weight of each particle is set as 1/L, and the weight is 
recursively updated based on the likelihood of the observation zi,t. To 
address the degeneracy problem, resampling described in [44] is applied 
in each step to eliminate the particles with small weights and duplicate 
the particles with high weights. After that, we can obtain a new set of 
particles y(p)i,t+1, p = 1, 2, …, L, with the same weight. The parameters of 
GMM can then be updated by a one-step EM algorithm [45] 
γ(p)i,t+1,k =
βi,t,kg
(
y(p)i,t+1|μi,t,k ,Σi,t,k
)
∑K
k=1βi,t,kg
(
y(p)i,t+1|μi,t,k ,Σi,t,k
), (8) 
μi,t+1,k =
(
∑L
p=1
γ(p)i,t+1,k
)− 1
∑L
p=1
y(p)i,t+1γ
(p)
i,t+1,k , (9) 
Σi,t+1,k =
(
∑L
p=1
γ(p)i,t+1,k
)− 1
∑L
p=1γ(p)i,t+1,k
(
y(p)i,t+1 − μi,t+1,k
)(
y(p)i,t+1 − μi,t+1,k
)T
,
(10) 
βi,t+1,k =
1
L
∑L
p=1
γ(p)i,t+1,k. (11) 
When the parameters of GMM at the time step (t+1) are available, 
the residual signal at the time step (t+1) can be defined with the updated 
parameters of Gaussian components as [46] 
ri(t+ 1) =
∑K
k=1
P
(
Ci,k
⃒
⃒zi,t+1
)
PMn
(
zi,t+1
)
, 0 < ri(t+ 1) < 1, (12) 
with 
P
(
Ci,k
⃒
⃒zi,t+1
)
=
βi,t+1,kg
(
zi,t+1
⃒
⃒λi,t+1,k
)
∑K
k=1βi,t+1,kg
(
zi,t+1
⃒
⃒λi,t+1,k
), (13) 
PMn
(
zi,t+1
)
= P
{
Dr
(
yi,t+1,Ci,k
)
≤ Dr
(
zi,t+1,Ci,k
)}
, (14) 
where P(Ci,k|zi,t+1) denotes the posterior probability of the monitored 
zi,t+1 belonging to the k-th Gaussian component, and PMn(zi,t+1) is the 
Mahalanobis distance of zi,t+1 from the center μi,t+1,k, which can be 
estimated as 
Dr
(
zi,t+1,Ci,k
)
=
(
zi,t+1 − μi,t+1,k
)T ( Σi,t+1,k + εI
)− 1( zi,t+1 − μi,t+1,k
)
, (15) 
where ε is a small positive number to avoid an ill-conditioned covariance 
matrix. 
Since Dr follows χ2 distribution with an appropriate degree of 
freedom, the threshold of the monitored residual signal should be set as 
the value of the confidence level (1-α), then the probability of monitored 
signals exceeding the alarm threshold for the i-th process node at time 
step (t+1) can be defined as 
PMONi (t+ 1) = φ
(
ri(t + 1) − 1 + α
σm
)
=
∫ ri(t+1)
− ∞
1
̅̅̅̅̅
2π
√
σm
e
(r− 1+α)2
2σ2m dr, (16) 
where σm = (1-α)/3 according to the three-sigma rule. It is worth noting 
that the monitored residuals larger than the threshold (1-α) should be 
scaled to ensure that when ri(t+1) is equal to 1, PMONi (t + 1)is also 
approximately equal to 1. The whole procedure for estimating the 
probability of process data exceeding the alarm threshold is summarized 
in Algorithm 1 below. 
3.2. Integrated real-time risk analysis 
When the temporal dependencies are introduced, the probability 
function of the top event (accident) can be transformed to the dynamic 
form from Eq. (1), and the accident probability at time tk is predicted as 
PTE(tk) = fFT(PBE1 (tk),PBE2 (tk),…,PBEN (tk)). (17) 
As shown in Fig. 4, DBN is used to present the dynamic behavior of 
risk models. Each BN node has two functional states, i.e., Occurrence 
(Occ) and Non-occurrence (Non). Since the state of basic events is 
observed firmly based on data monitoring at each time, the conditional 
probability distribution of P(BEi(tk)|BEi(tk − 1)) is neglected [47]. 
Therefore, the risk index (posterior probability) of each basic event can 
be estimated based on the Bayes rule 
P(BEi(tk) = Occ|TE = Occ)
=
P(TE = Occ|BEi(tk) = Occ)PBEi (tk)∑
BEi(tk)=Occ,Non
P(BEi(tk))P(TE = Occ|BEi(tk))
, i = 1, 2,…N. (18) 
From the above equations, we can see that the key point of PRA is to 
estimate the probability of BEi at time tk. In this paper, a BN (as shown in 
the red frame in Fig. 4) is developed to describe the dependencies be-
tween prior probability PBEi and the corresponding alarm probability 
PMONi for the i-th event, and then we can haveP(BEi, MONi) =
P(MONi|BEi)PBEi . Since PMONi obtained in Section 3.1 is a continuous 
value between 0 and 1, the continuous probability needs to be dis-
cretized to construct the CPT of node MONi. Assume that the state of the 
node MONi is divided into Q discretization states S = {S1, S2, …, SQ}, 
where Q should be an even number and the state SQ corresponds to the 
probability from (Q − 1)/Q to 1. According to the definition in Section 
3.1, the probability PMON is 0.5 for a given point on the alarm threshold. 
As shown in Fig. 4, we assume that the monitored signal of an event 
obeys a normal distribution, in which its probability density function 
(PDF) is denoted byf(⋅); the upper and lower limits of thef(⋅) are denoted 
by TU and TL, respectively; and τ is the alarm threshold set for the 
Algorithm 1 
Probability estimation for process nodes based on GMM and PF 
Offline 
Inputs: historical fault-free data yi 
Outputs: Λi = {βi, μi, Σi} 
Step: Train GMM using the EM algorithm. 
Online (at time t+1) 
Inputs: Λi,t = {βi,t, μi,t, Σi,t}, zi,t, zi,t+1 
Outputs: Λi,t+1 = {βi,t+1, μi,t+1, Σi,t+1}, PMONi(t+1) 
Step 1: Generate a set of random particlesy(p)i,t , p = 1, 2, …, L, according to the Λi,t =
{βi,t, μi,t, Σi,t}. 
Step 2: Calculate the weight of particles based on zi,t. 
Step 3: Resample to eliminate particles with small weights and obtain a new set of 
particlesy(p)i,t+1, p = 1, 2, …, L. 
Step 4: Update the GMM parameters by Eq. (8)-(11) and get Λi,t+1 = {βi,t+1, μi,t+1, 
Σi,t+1}. 
Step 5: Calculate the monitored residual ri(t+1) based on Λi,t+1 and zi,t+1 by Eq. (12)- 
(15). 
Step 6: Calculate the probability PMONi(t+1) by Eq. (16) 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
6
monitored signal. If the state of basic event is Non, the monitored signal 
will be in the range from TL to τ (corresponding to PMON from 0 to 0.5) 
with the probability of 
∫ τ
TL
f(⋅) and will be in the range from τ to TU 
(corresponding to PMON from 0.5 to 1) with the probability of 
∫ TU
τ f(⋅), 
and vice versa. Using this idea, the distribution of the monitored signal 
can be scaled into CPT of BN. The conditional probability P(Sj |BE) can 
be estimated with a pre-set threshold τ and PDFf(⋅)as 
P
(
Sj
⃒
⃒BE
)
=
⎧
⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎩
∫ 2τj/Q
2τ(j− 1)/Q
f (x)dx
∫ TU
TL
f (x)dx
, 0 <j ≤
Q
2
∫ τ+(2j− Q)(TU − τ)/Q
τ+(TU − τ)(2j− Q− 2)/Q
f (x)dx
∫ TU
TL
f (x)dx
,
Q
2
< j ≤ Q
. (19) 
When the historical data is available, for example, when the histor-
ical fault-free data is useful, the kernel density estimation (KDE) can be 
used to estimate the PDF of the residual signal under normal conditions, 
and the estimated PDF can be applied as the f(⋅) for P(Sj
⃒
⃒BE = Non). The 
∑j=Q
j=Q/2+1P(Sj
⃒
⃒BE= Non) can denote the false alarm rate (FAR), and the 
∑j=Q
j=Q/2+1P(Sj
⃒
⃒BE= Occ) can represent the fault detection rate (FDR). 
When no historical data is available, thef(⋅)can be defined as the PDF of 
normal distribution N(μ,σ), and the TU and TL can be set as (μ+3σ) and 
(μ-3σ), separately. The value of μ and σ can be adjusted arbitrarily to 
ensure TL = 0. A high value of FAR indicates many false alarms will be 
triggered when the process is normal, while a low FDR implies few 
alarms can be raised when the process is abnormal. The threshold τ is set 
to determine the value of FAR and FDR, which can represent the un-
certainty in the risk monitoring model. For the case of BE = Non, the 
threshold τ can be set as (μ+2.5σ) to guarantee nearly 95% alarms are 
correct for the fault-free case; for the case of BE = Occ, the threshold τ 
can be set as (μ-2.5σ) to ensure only 5% alarms can not be triggered for 
the faulty case. When PMONi is observed, the corresponding state of the 
node MONi, denoted by S∗i , can be obtained, and PBEi (tk) can then be 
computed as follows 
PBEi (tk) =
P
(
S∗i
⃒
⃒BEi
)
PBEi
∑
BEi=Occ,NonP(BEi)P
(
S∗i
⃒
⃒BEi
). (20) 
The CPT in Fig. 4 shows an example where Q and μ are 10 and 3, 
respectively. As an example, for a basic event, if the alarm probability 
PMON is 0.25 and prior probability PBE is 0.05, the S∗ will be S3, and the 
updated probability will be 9.14× 10− 4 × 0.05/(9.14 × 10− 4 × 0.05 +
0.4072 × 0.95) = 1.18× 10− 4. If PMON is less than 0.5, the updated 
probability will be less than the prior probability. 
Furthermore, based on the idea of discretization, the P(MONi(tk)| 
BEi) at time tk can be computed directly by 
P(MONi(tk)|BEi) =
⎧
⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎩
∫ 2τPMONi (tk)+Δ
2τPMONi (tk)− Δ
f (x)dx
∫ TU
TL
f (x)dx
, 0 ≤ PMONi (tk) ≤ 0.5
∫ 2τ− TU+2PMONi (tk)(TU − τ)+Δ
2τ− TU+2PMONi (tk)(TU − τ)− Δ
f (x)dx
∫ TU
TL
f (x)dx
, 0.5 < PMONi (tk) ≤ 1
,
(21) 
where Δ is a small positivenumber representing the discrete interval. By 
Eq. (21), we can directly input the value of PMONi into the model without 
identifying its corresponding state S∗i , so that we can get continuously 
updated probability based on different PMONi . It is worth mentioning that 
both the first and the second equations of P(MONi(tk)|BEi) will be equal 
to zero when PMONi (tk) is equal to zero or one. Therefore, in this work, we 
compress the boundary of PMONi (tk) to a new boundary from 0.2 to 0.8 by 
PMONi (tk) = 0.6× PMONi (tk)+ 0.2. The figure on the left side of CPT in 
Fig. 4 (as shown in the blue frame) shows the value of P(MON(tk)| BE) 
with Δ equal to 0.2. Finally, thePBEi (tk)can be computed based on the 
value ofP(MONi(tk)|BEi = Occ)andP(MONi(tk)|BEi = Non) as follows 
Fig. 4. Diagram of the proposed PRA method. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
7
PBEi (tk) =
P(MONi(t)|BEi = Occ)PBEi∑
BEi=Occ,NonP(BEi)P(MONi(t)|BEi)
. (22) 
4. Implementation and model discussion 
In this section, two case studies are used to demonstrate the benefits 
of the proposed method brought to PRA. In Section 4.1, an experimental 
study with comparative analysis is presented. Then in Section 4.2, data 
monitoring-based PRA is performed for MPD operations to evaluate the 
proposed model. 
4.1. Case study 1: RT 580 experimental setup 
4.1.1. System description 
In the first case, the RT 580 experimental setup (shown in Fig. 5) 
[48], situated at the Memorial University, is applied to demonstrate the 
proposed method. The RT 580 setup controls the level of process tank 
(B2) at 40% through a circuit with a collecting tank (B1), pump (P1), 
and valves. The pneumatic control valve (V7) is used as the actuator. 
The FT shown in Fig. 6 represents the accident mechanism of “dry-out in 
tank B2”. Eight basic events (BE1 ~ BE8) are considered, and their prior 
probabilities can be found in [48]. 
The flow rate and level in the process tank are monitored during the 
experiment. The basic event BE1 “Sensor is broken” is introduced into 
the system. Fig. 7 shows the experimental data collected from the stable 
operation to the failure of the system. When the level sensor is broken, 
the controller acquires the default value (high value) of the level in the 
tank. The controller attempts to close the control valve (V7), causing the 
flow and level in tank B2 to drop. Since the system operates in a signal 
condition, we consider the node BE1 as an environmental node. Ac-
cording to Ref. [48], a flow rate less than 100 L/h and a level less than 
10% will cause a dry-out situation. 
Therefore, the probability of BE1 is determined by flow rate and level 
monitoring, and the normal reference value and the threshold value of 
flow rate and level are set as {400, 40} and {100, 10}, respectively. 
4.1.2. PRA results and discussions 
The figure on the left-hand side of Fig. 8 shows the probability of BE1 
signals exceeding alarm thresholds. The real-time probability of BE1 
shown in the middle figure of Fig. 8 can be calculated by integrating its 
prior probability and alarm probability. The first two figures in Fig. 8 
present probability profiles with similar growth trends. As defined in 
Section 3.1, a flow rate less than 100 L/h and a level less than 10% will 
result in the alarm probability of BE1 higher than 0.5, indicating that the 
abnormality of its monitored signal keeps growing and may lead to 
abnormal events. Although the alarm probability by the end of BE1 is as 
high as a probability close to 0.9, the real-time probability of BE1 is still 
small due to the low prior probability. In this case, the real-time prob-
ability of BE1 is consistent with its prior risk knowledge. The figure on 
the right-hand side of Fig. 8 shows the real-time probability of dry-out 
event calculated by the proposed method and dynamic fault tree 
(DFT) [48]. In DFT, the probability of BE1 (sensor is broken) is equal to 1 
when the flow rate is less than 100 L/h or the level is less than 10%, 
which will result in a real-time probability of 1 for the dry-out event. 
However, when the process variable exceeds the alarm threshold, the 
corresponding accident may not always occur. Alarms can also be 
caused by process disturbances or signal noise. In contrast, the accident 
probability predicted by the proposed method does not increase 
immediately when the process variables exceed the alarm threshold. The 
probability curves of the dry-out accident undertaken with the proposed 
method are smoother and show significant variations in values 
compared to the normal case, and thus are more reasonable. 
To search for critical basic events, the normalized ratios of posterior 
probabilities are compared. Fig. 9 presents the static posterior estimates 
of basic events in the RT 580 setup when the dry-out event occurs. 
Before the level sensor is damaged, static posterior analysis can find 
some basic events with potential dangers of the dry-out accident, such as 
BE1, BE7, and BE8. Fig. 10 compares posterior estimates of basic events in 
the event of a dry-out accident between the proposed method and DFT. 
Since BE1 has been monitored and the monitoring results show that BE1 
is abnormal when the accident occurs, the posterior probability ratio of 
BE1 increases rapidly and is significantly higher than other basic events. 
Both of the proposed method and DFT can find the most critical event in 
real-time. However, on the basis of posterior results by DFT, the pro-
posed method can further find that the posterior probability of BE8 is 
slightly higher than the other six events, indicating that BE8 still needs 
further attention. 
4.2. Case study 2: managed pressure drilling operations 
4.2.1. System description 
Managed pressured drilling (MPD) system is a typical closed-loop 
kick control system used in deepwater drilling operations. As shown in 
Fig. 11, MPD utilizes the choke and the pump to adjust annular pressure 
to meet the requirements of bottomhole pressure (BHP) and observes the 
BHP condition based on the measurements of flow and hydrodynamic 
pressure. Fig. 12 shows the FTA of MPD operations. The top event (kick) 
results from the combined effect of both the abnormal BHP and the 
failure of the MPD system [50]. Twenty-one basic events (BE1 ~ BE21) 
are listed with the potential accidents. The prior probability PBEi, i = 1, 2, 
…, 21 of basic events can be found in [50]. 
In this case, all basic events are divided into two groups; the first 
group, 
GBE = {BE7,BE8,BE12,BE14,BE16,BE18}, (23) 
represents the events that can be updated in real-time, in which the 
probability of BE12 and BE14 can be set as 1 when the corresponding 
signal is lost, whereas the probabilities of the remaining basic events are 
assumed to be constant values equal to their prior probabilities. In 
addition, the BE7 (abnormal pressured zone) and BE8 (loss of circulation) 
are regarded as environmental nodes, and their probabilities are upda-
ted based on fluid and pressure monitoring. The BE16 (rig pump failure) 
and BE18 (backpressure pump failure) are considered process nodes 
since their reference signal values and alarm thresholds can be affected 
by changes in pump operating conditions. 
Data collected from MPD in conjunction with experimental equip-
ment data are used in the second case study. The fluid and pressure Fig. 5. The simplified flow diagram of RT 580 experimental setup. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 1087008
measurement data come from the actual drilling of Well GS19 in the 
Qixia Formation [51]. The density of injected and returned drilling fluid 
is measured with a flowmeter to identify overflow or lost circulation, 
Fig. 6. FT for “Dry-out in tank B2” incident (adapted from the Ref. [48]). 
Fig. 7. Flow rate and level measurements during the experiment [48]. 
Fig. 8. Forward PRA updating of BE1 and dry-out event. 
Fig. 9. Static posterior estimation of RT 580 setup. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
9
and the automatic choke manifolds are used to regulate wellhead back 
pressure to maintain the BHP in a safe range. Fig. 13 shows the MPD 
data, e.g., casing pressure, inlet density, outlet density, and BHP 
equivalent circulation density, collected from the drilling depth of 4035 
m to 4235 m. Two layers with lost circulation were encountered at a 
depth of 4145.13 m and 4232.22 m, while the lost layers were plugged 
successfully through the automatic injection of bridging mud at 
managed pressure. In addition, it can be seen from the change in 
pumping rate that the operating condition of the rig pump occurs at 
depths of 4140 m, 4210 m, and 4230 m. Observed from the dataset, 
higher casing pressure and BHP equivalent circulation density indicate a 
high risk of BE7, and the outlet density increases and is slighter higher 
than the inlet density when the loss of circulation occurs. In practice, a 
basic event can be modeled based on multiple signals using similar 
approaches. 
The health status of oil pumps can be represented by the vibration 
features of the bearing inside the pump, and thus, the vibration data of 
the bearing provided by Case Western Reserve University (CWRU) 
Bearing Data Center [52] is used to test the monitoring ability of the 
proposed method. The database includes normal bearing monitoring 
data in four conditions, labeled from Normal 1 to Normal 4, corre-
sponding to the motor speed from 1797 rpm to 1730 rpm. We assume 
that 10000 vibration data points can be obtained for every 10 m step of 
drilling depth to match the drilling process in Fig. 13. To simplify the 
Fig. 10. Posterior estimation of RT 580 setup at time 180 s. 
Fig. 11. MPD Process Flow Diagram (adapted from the Ref. [49]). 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
10
model, only the vibration data of the drive end acceleration are used. 
The rig pump is assumed to be fault-free during the whole drilling 
process but has multiple working conditions up to the actual case, 
whereas for the backpressure pump, an abrupt rolling element fault is 
assumed to occur at a depth of 4200 m. In this case, the rest of the 
Normal 1 and Normal 4 data are used to train the monitoring model for 
the rig pump and backpressure pump, respectively. The use of vibration 
data is shown in Fig. 14. More details of the CWRU database can be 
found in [53]. 
4.2.2. Real-time PRA results 
In this section, the PRA of MPD operations is performed following the 
procedure in Fig. 2. According to the actual drilling of Well GS19 [51], 
the BHP is kept at a slight over-balance during drilling, and the annular 
back pressure and BHP equivalent circulation density are controlled at 
2–5 MPa and 2.46–2.52 g/cm3. Therefore, the normal reference value 
and the threshold value of casing pressure, BHP equivalent circulation 
density, and the difference between outlet and inlet densities are set as 
{3.5, 2.49, -0.03} and {5, 2.52, 0}, respectively. In addition, three sta-
tistical time-domain features are extracted from the vibration 
measurements 
⎧
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨
⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎪⎩
x1 =
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
1
N
∑N
k=1
y2k
√
√
√
√
x2 =
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
1
N − 1
∑N
k=1
(yk − y)2
√
√
√
√
x3 =
̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
1
N
∑N
k=1
y2k
√
1
N
∑N
k=1
|yk|
, (24) 
where x1 is the root mean square, x2 is the standard deviation, and x3 is 
the shape factor [54]. In Eq. (24), y is the vibration signal, and N is the 
number of vibration points equal to the product of time interval length 
and sampling frequency. In this case, N is equal to 100. The EM algo-
rithm is used to train the GMM based on the historical vibration data, 
and then we can obtain 
Fig. 12. FTA for kick control incident (adapted from the Ref. [50]). 
βBE16 = [ 0.038 0.070 0.264 0.629 ],
μBE16 =
⎡
⎢
⎣
0.064 0.087 0.089 0.065
0.057 0.084 0.088 0.063
1.266 1.221 1.179 1.230
⎤
⎥
⎦,
Σ(1)BE16 =
⎡
⎢
⎣
5.64e − 5 6.24e − 5 5.33e − 5
6.24e − 5 8.31e − 5 1.50e− 4
5.33e − 5 1.50e− 4 2.30e− 3
⎤
⎥
⎦,Σ(2)BE16 =
⎡
⎢
⎣
1.20e − 4 1.19e − 4 − 5.95e − 6
1.19e − 4 1.27e − 4 2.21e− 5
− 5.95e − 6 2.21e− 5 7.60e − 4
⎤
⎥
⎦
Σ(3)BE16 =
⎡
⎢
⎣
9.28e − 5 9.04e − 5 − 1.27e − 5
9.04e − 5 9.09e − 5 − 1.47e − 5
− 1.27e − 5 − 1.47e − 5 4.49e − 4
⎤
⎥
⎦,Σ(4)BE16 =
⎡
⎢
⎣
1.22e− 4 1.21e − 4 − 7.83e − 5
1.21e − 4 1.22e − 4 − 8.20e − 5
− 7.83e − 5 − 8.20e − 5 1.40e − 3
⎤
⎥
⎦,
(25) 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
11
where ΛBE16 = {βBE16, μBE16, ΣBE16} and ΛBE18 = {βBE18, μBE18, ΣBE18} are 
the trained GMM for BE16 and BE18, respectively. 
After determining model parameters, the probability PMON7 and 
PMON8 that online signals exceed alarm thresholds are calculated for BE7 
and BE8 based on fluid and pressure monitoring by Section 3.1.1. In 
addition, with the assumption that the interval length is equal to 1 m, the 
alarm probabilities PMON16 and PMON18 are computed in real-time for BE16 
and BE18 based on the trained GMMs and the dynamic updating method 
shown in Algorithm 1. Fig. 15 exhibits the results of real-time proba-
bility PMON for four basic events. As depicted in Fig. 15 (a), when the lost 
circulation occurs at a depth of 4145.13 m and 4232.22 m, PMON8 (Loss 
of circulation) increases significantly. However, since this case only 
applies a single variable to evaluate the lost circulation, the measure-
ment noise will seriously affect the evaluation results, so that sometimes 
under normal working conditions, BE8 also has a higher alarm proba-
bility PMON8 . In contrast, PMON7 (abnormal pressured zone) is assessed by 
two variables and thus has higher stability. As shown in Fig. 15 (b), 
PMON16 (rig pump failure) will not change significantly due to the 
changes in operating conditions. Meanwhile, the proposed method can 
effectively monitor the fault of BE18 (backpressure pump failure) at a 
depth of 4200 m. 
Next, the estimated alarm probabilities PMONi are integrated with 
prior probabilities to evaluate the risk of dynamic nodes in real-time. 
According to Eq. (21) in Section 3.2, the conditional probability P 
(MONi| BEi) at each depth can be computed based on PMONi , and P(MONi| 
BEi) is incorporated into prior probabilities to calculate the real-time 
risks of dynamic nodes by Eq. (22). The obtained real-time probability 
of basic events BE7, BE8, BE16, and BE18 is shown in Fig. 16. In this case, 
the model parameter Δ is set to 0.2, and μ is set to 3. Detailed parameter 
settings will be discussed further in Section 4.2.4. It can be observed that 
if the alarm probability PMON obtain in Fig. 15 is less than 0.5, the risk of 
basic events in Fig. 16 will not change significantly. Fig. 16 shows that 
therisk of BE8 increases when lost circulation occurs at depths of 
Fig. 14. Vibration measurements for the drilling process. 
Fig. 13. Measurements during the MPD of Well GS19 [51]. The BHP equivalent 
circulation density is only monitored at depths of 4090 m and 4170 m. 
βBE18 = [ 0.046 0.074 0.053 0.162 ],
μBE18 =
⎡
⎢
⎣
0.062 0.062 0.075 0.074
0.055 0.061 0.073 0.070
1.238 1.231 1.290 1.207
⎤
⎥
⎦,
Σ(1)BE18 =
⎡
⎢
⎣
4.93e − 5 3.79e − 5 2.60e − 5
3.79e − 5 3.83e − 5 2.88e− 5
2.60e − 5 2.88e− 5 1.10e− 3
⎤
⎥
⎦,Σ(2)BE18 =
⎡
⎢
⎣
7.97e − 5 7.81e − 5 − 1.89e − 5
7.81e − 5 7.93e − 5 − 1.67e− 5
− 1.89e − 5 − 1.67e− 5 9.94e − 4
⎤
⎥
⎦
Σ(3)BE18 =
⎡
⎢
⎣
1.28e − 4 1.24e − 4 − 1.99e − 5
1.24e − 4 1.25e − 4 − 1.94e − 5
− 1.99e − 5 − 1.94e − 5 1.20e− 3
⎤
⎥
⎦,Σ(4)BE18 =
⎡
⎢
⎣
7.99e− 5 8.22e − 5 − 3.64e − 5
8.22e − 5 9.47e − 5 − 3.23e − 5
− 3.64e − 5 − 3.23e − 5 1.00e − 3
⎤
⎥
⎦,
(26) 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
12
4145.13 m and 4232.22 m, and the risk of BE18 rises when pump fault 
occurs at a depth of 4200 m. Therefore, the proposed method is vali-
dated to be effective. 
When we get the risk of all basic events at a certain drilling depth, the 
probability of kick can be estimated by Eq. (17). As shown in Fig. 17 (a), 
a single lost circulation accident will not significantly increase the 
probability of a kick due to the role of the MPD system as a safety barrier. 
When the MPD system can work normally, the closed-loop control of 
MPD can balance the BHP to alleviate the lost circulation and avoid the 
kick, and also, this phenomenon is entirely consistent with the actual 
drilling case [51]. In contrast, the MPD pump equipment failed at a 
depth of 4200m and could not adjust the BHP normally, and thus the 
probability of kick increased instantly. At a depth of 4232 m, the lost 
circulation and the backpressure pump fault occur at the same time, and 
therefore, the probability of kick reaches the highest value, which also 
indicates that the kick is a result of the combined effect of both the 
abnormal BHP and the failure of the MPD system. As a comparison, 
Fig. 17 (b) shows the kick probability calculated without considering the 
prior probability of basic events, which is a similar idea that has been 
widely used in many reported works [33–37]. Although the comparative 
Fig. 16. The real-time risk of basic events. 
Fig. 15. The probability of online signals exceeding alarm thresholds. 
Fig. 17. The real-time probability of kick. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
13
method can successfully identify the rise in risk after the depth of 4220 
m, the obtained kick risks have large fluctuations due to measurement 
noise, even when the MPD system is working normally (from a depth of 
4035 m to 4150 m). 
Beyond the real-time probability analysis of kick, it is also essential 
to discover the critical events with potential kick risk. The posterior risk 
estimation can be performed based on obtained real-time probabilities 
by Eq. (18) and the proposed method. In this case, the kick is adopted as 
the top event to calculate the posterior probability of basic events. 
Fig. 18 presents the normalized posterior estimates at some 
representative drilling depths. The first graph in Fig. 18 shows the result 
of the static posterior analysis. Before drilling operations, static poste-
rior analysis can find some basic events with potential dangers of kick 
accident, such as BE1, BE2, BE4, BE7, BE15, and BE16. However, since BE7 
and BE16 have been monitored and the monitoring results show that 
they are normal, the posterior probability ratio of the events BE7 and 
BE16 decreases at a depth of 4060 m. At a depth of 4150 m, the lost 
circulation occurs, and thus, the BE8 (loss of circulation) increases, 
indicating that the BE8 is quite critical and should be noted at this depth. 
From the depth of 4150 m to 4210 m, the lost layers are successfully 
Fig. 18. Real-time posterior probability of MPD operations using the proposed method. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
14
plugged, but at a depth of 4210 m, the backpressure pump malfunction 
occurs and the posterior probability ratio of the events BE18 increases 
correspondingly. In addition, with the rise of the posterior probability of 
BE18, the posterior probability of BE19 also increases, and thus BE19 and 
BE18 may have similar accident propagation paths. A similar analysis 
can be conducted for the last two graphs. At a depth of 4230 m, the rig 
pump changes its operating conditions to suppress the lost circulation so 
that the posterior probability of BE8 and BE16 is slightly increased. At a 
depth of 4235 m, the rig pump returns to its initial operating condition, 
and the corresponding posterior probability drops. For comparison, 
Fig. 19. Real-time posterior probability of MPD operations without considering the prior probability. 
Fig. 20. The probability of BE16 signals exceeding alarm thresholds. 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
15
when prior risk is not considered, the posterior risks of BE8 (loss of 
circulation) are still high, as shown in Fig. 19, even if the layers are not 
lost, or the lost layers are successfully plugged (at depths of 4060 m, 
4230 m, and 4235 m). In addition, at a depth of 4150 m, the comparative 
method incorrectly gets a high posterior risk of BE16 since the prior risk 
of rig pump failure is not incorporated in the estimation. Overview, the 
static risk analysis is hard to capture the deviations and changes in MPD 
operations, and the risk cannot be well estimated if only the monitoring 
information is considered without considering prior risks. However, 
once the prior risks are to be considered, CPTs need to be set up to 
connect prior probabilities to data monitoring, which often requires 
historical accident data or expert knowledge in reported works. There-
fore, the proposed real-time PRA is particularly effective in evaluating 
the time-dependent behavior of MPD operations in the presence of data 
scarcity and missing expertise. 
4.2.3. The comparison of GMM and GMM-PF 
MPD is a typical intermittent operation system. The operating con-
dition of the rig pump and backpressure pump will change over time, 
and the transient behavior can occur during phase shifts. GMM can be 
utilized to construct the monitoring model for process nodes. However, 
as the condition of mechanical equipment may change in MPD opera-
tions, the component parameters of GMM require updates over time. 
Therefore, PF is developed to update the GMM in Algorithm 1. This part 
is applied to demonstrate the positive impact of particle filtering-based 
parameter updating. 
As shown in Fig. 14 (a), the operating condition of the rig pump will 
change during the drilling. Since the features in Eq. (24) cannot distin-
guish different operating condition, we build a simple neural network 
model to extract a new set of three featuresfrom vibration signals, which 
can effectively present the difference between different conditions, and 
the input of the neural network is 100 vibration points. New features are 
used to train the GMM for the rig pump (BE16) and the backpressure 
pump (BE18). The probability of BE16 and BE18 signals exceeding alarm 
thresholds based on GMM are shown in Fig. 20 (a) and Fig. 21 (a), 
respectively. When we use PF to update the GMM parameters, the ob-
tained alarm probability of BE16 and BE18 are shown in Fig. 20 (b) and 
Fig. 21 (b). If the alarm probability exceeds 0.5, it means that the cor-
responding mechanical equipment may be malfunctioning. As presented 
in Fig. 20 (a), the GMM-based monitoring model will incorrectly eval-
uate other normal operating conditions as faults. Specifically, the motor 
speed of the rig pump changes from 1797 rpm to 1750 rpm at a depth of 
4140 m and increases from 1750 rpm to 1772 rpm at a depth of 4210 m. 
Fig. 21. The probability of BE18 signals exceeding alarm thresholds. 
Fig. 22. Sensitivity analysis of μ 
Fig. 23. Sensitivity analysis of Δ 
R. He et al. 
Reliability Engineering and System Safety 226 (2022) 108700
16
Although the rig pump still works normally, the GMM-based monitoring 
model incorrectly identifies the change in operating conditions as a fault 
since most of the corresponding probabilities exceed 0.5 after depths of 
4140 m and 4210 m. In contrast, as shown in Fig. 20 (b), PF can effec-
tively adjust the GMM parameters to a reasonable range so that GMM-PF 
will not be too sensitive to changes in pump operating conditions. In 
addition, at a depth of 4200 m, the backpressure pump experienced a 
sudden rolling element fault, as shown in Fig. 14 (b). As presented in 
Fig. 21, both GMM and GMM-PF can identify the fault of the back-
pressure pump at a depth of 4200 m. Overall, the proposed GMM-PF 
method is more applicable when the monitoring parameters are 
affected by varying working conditions. 
4.2.4. Model parameter analysis 
This work contains three models presented in Section 3. The first 
monitoring model shown in Section 3.1.1 is a deterministic model, 
which only requires reference values and thresholds for process mea-
surements instead of any pre-set parameters. For the GMM-PF model 
presented in Section 3.1.2, the parameter settings of the model can refer 
to many existing references, such as Ref. [42] and [43]. As a result, this 
section focuses on discussing the parameter setting of the integrated PRA 
model shown in Section 3.2. 
We need to set only two parameters, μ and Δ, for the proposed PRA 
model, and other parameters can be derived from these two parameters. 
This section provides a sensitivity analysis to analyze how to determine 
model parameters for PRA. First, we set Δ to 0.05 and prior probability 
to 0.05 to observe the influence of the change of parameter μ on the 
posterior probability. As shown in Fig. 22, whether the PMON value is 0.4 
or 0.6, the posterior probability will remain stable when μ is greater than 
about 2, and thus we can set the parameter μ to a number greater than 2. 
Fig. 23 shows the sensitivity analysis result of parameter Δ, in which the 
μ is set as 3 and the prior probability is set as 0.05. Δ is a small number 
denoting the discrete interval. When the parameter Δ is less than 0.3, the 
posterior probabilities are stable, and therefore, we can set Δ to a very 
small positive number for the PRA model. 
4.2.5. Model uncertainty and validation 
The uncertainty and validation of risk analysis models have been 
widely studied. The model uncertainty may greatly impact the estima-
tion accuracy of the probability estimation of basic events, and then 
affect the calculated risk. One way to represent the uncertainty in PRA is 
to conduct a bounding analysis by using probability boxes (P-boxes) 
[55]. In the bounding analysis, the prior probability of basic events can 
be assessed as a credible interval rather than a single point. By propa-
gating the uncertainty, a credibility interval can also be obtained for the 
estimated kick risk and posterior risks, which can help decision-makers 
understand the confidence of risk estimations. 
Three-axiom-based validation [56] can be performed to give a partial 
validation for the proposed PRA model. If the proposed model is robust, 
the estimated probability of basic events is often sensitive but rarely 
presents abrupt variations due to any slight changes in corresponding 
alarm probability. As shown in Fig. 20 and Fig. 21, if PMON is less than 
0.5, the real-time risk of the corresponding basic event will only change 
slightly. On the contrary, when PMON is greater than 0.5, the real-time 
risk of the corresponding basic event will increase significantly. With 
a minor change in the alarm probability, the real-time risk of basic 
events is sensitive but rarely exhibits abrupt changes. In addition, once 
the monitoring results show that the event status is abnormal, the 
real-time risk of the basic event will rise reasonably. The sensitive 
analysis results satisfy the axioms and thus partially validate the pro-
posed PRA model. 
Conclusions 
This paper proposes a novel integrated PRA method for petrochem-
ical industries based on data monitoring. The proposed method can fully 
satisfy the timing requirement of process operations. Additionally, the 
proposed probability estimation not only considers the update of the 
failure probability of safety barriers but also assesses the probability of 
abnormal events. Furthermore, PF-based GMM can estimate the alarm 
probability of process signals effectively during phase shafts in industrial 
systems. By comparing different risk analysis methods, the main con-
clusions from the case study validate the effectiveness and necessity of 
the proposed PRA method. 
The significant contribution of the proposed approach is the inte-
gration of data monitoring with DBN for PRA in real-time. Specifically, 
the probability of the monitored signal exceeding the alarm threshold is 
estimated first for environmental and process monitoring based on the 
multivariate method and PF-based GMM, respectively. Data monitoring 
can be effectively incorporated into PRA in real-time by the proposed 
DBN with the estimated alarm probability. The proposed DBN model 
will not require additional expert knowledge and historical accident 
data to determine the conditional relationship between data monitoring 
and prior risk so that is more suitable for practical engineering appli-
cations. We plan to collect actual process data to test the model and 
incorporate effective monitoring models into PRA. 
CRediT authorship contribution statement 
Rui He: Conceptualization, Methodology, Writing – original draft, 
Writing – review & editing. Jingyu Zhu: Investigation, Writing – review 
& editing. Guoming Chen: Writing – review & editing, Funding 
acquisition. Zhigang Tian: Writing – review & editing. 
Declaration of Competing Interest 
No conflict of interest exists in the submission of this manuscript, and 
all authors approve the manuscript for publication. We would like to 
declare on behalf of co-authors that the work described was original 
research that has not been published previously, and not under 
consideration for publication elsewhere, in whole or in part. We declare 
that we do not have any commercial or associative interest that repre-
sents a conflict of interest in connection with the work submitted. 
Acknowledgments 
The authors greatly appreciate the help from the Centre for Offshore 
Engineering and Safety Technology (COEST) at China University of 
Petroleum (East China), and gratefully acknowledge the financial sup-
port provided by China National Key Research and Development Pro-gram (No: 2016YFC0304005), and National Natural Science Foundation 
of China (No: 52004142). 
References 
[1] Khan F, Abbasi S. Techniques and methodologies for risk analysis in chemical 
process industries. J Loss Prevent Proc 1998;11(4):261–77. 
[2] Hu Y, Parhizkar T, Mosleh A. Guided simulation for dynamic probabilistic risk 
assessment of complex systems: Concept, method, and application. Reliab Eng Syst 
Safe 2022;217:108047. 
[3] Khan F, Abbasi SA. Analytical simulation and PROFAT II: a new methodology and a 
computer automated tool for fault tree analysis in chemical process industries. 
J Hazard Mater 2000;75(1):1–27. 
[4] Ferdous R, Khan F, Sadiq R, Amyotte P, Veitch B. Handling data uncertainties in 
event tree analysis. Process Saf Environ 2009;87(5):283–92. 
[5] Aqlan F, Mustafa Ali E. Integrating lean principles and fuzzy bow-tie analysis for 
risk assessment in chemical industry. J Loss Prevent Proc 2014;29:39–48. 
[6] Khakzad N, Khan F, Amyotte P. Dynamic safety analysis of process systems by 
mapping bow-tie into Bayesian network. Process Saf Environ 2013;91(1–2):46–53. 
[7] Moradi R, Groth KM. Modernizing risk assessment: A systematic integration of PRA 
and PHM techniques. Reliab Eng Syst Safe 2020;204:107194. 
[8] OREDA-offshore reliability data handbook. 4th ed. Trondheim: SINTEF; 2002. 
[9] Li X, Zhang Y, Abbassi R, Khan F, Chen G. Probabilistic fatigue failure assessment 
of free spanning subsea pipeline using dynamic Bayesian network. Ocean Eng 
2021;234:109323. 
[10] He R, Chen G, Wang Y, Jiang S, Zhi C. A quantitative risk analysis model 
considering uncertain information. Process Saf Environ 2018;118:361–70. 
R. He et al. 
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0001
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0001
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0002
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0002
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0002
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0003
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0003
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0003
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0004
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0004
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0005
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0005
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0006
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0006
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0007
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0007
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0008
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0009
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0009
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0009
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0010
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0010
Reliability Engineering and System Safety 226 (2022) 108700
17
[11] Liu Z, Ma Q, Cai B, Shi X, Zheng C, Liu Y. Risk coupling analysis of subsea blowout 
accidents based on dynamic Bayesian network and NK model. Reliab Eng Syst Safe 
2022;218:108160. 
[12] Meel A, Seider W D. Plant-specific dynamic failure assessment using Bayesian 
theory. Chem Eng Sci 2006;61:7036–56. 
[13] Kalantarnia M, Khan F, Hawboldt K. Dynamic risk assessment using failure 
assessment and Bayesian theory. J Loss Prevent Proc 2009;22:600–6. 
[14] Li X, Chen G, Khan F, Xu C. Dynamic risk assessment of subsea pipelines leak using 
precursor data. Ocean Eng 2019;178:156–69. 
[15] Yu H, Khan F, Veitch B. A flexible hierarchical Bayesian modeling technique for 
risk analysis of major accidents. Risk Anal 2017;37(9):1668–82. 
[16] Meng X, Li X, Wang W, Song G, Chen G, Zhu J. A novel methodology to analyze 
accident path in deepwater drilling operation considering uncertain information. 
Reliab Eng Syst Safe 2021;205:107255. 
[17] Yang X, Haugen S, Paltrinieri N. Clarifying the concept of operational risk 
assessment in the oil and gas industry. Safety Sci 2018;108:259–68. 
[18] Khakzad N, Khan F, Amyotte P. Quantitative risk analysis of offshore drilling 
operations: A Bayesian approach. Safety Sci 2013;57:108–17. 
[19] Pasman H, Reniers G. Past, present and future of Quantitative Risk Assessment 
(QRA) and the incentive it obtained from Land-Use Planning (LUP). J Loss Prevent 
Proc 2014;28:2–9. 
[20] Shalev DM, Tiran J. Condition-based fault tree analysis (CBFTA): A new method for 
improved fault tree analysis (FTA), reliability and safety calculations. Reliab Eng 
Syst Safe 2007;92(9):1231–41. 
[21] Liu Q, Dong M, Peng Y. A novel method for online health prognosis of equipment 
based on hidden semi-Markov model using sequential Monte Carlo methods. Mech 
Syst Signal Pr 2012;32:331–48. 
[22] Zeng Z, Zio E. Dynamic risk assessment based on statistical failure data and 
condition-monitoring degradation data. IEEE T Reliab 2018;67(2):609–22. 
[23] Gomes JPP, Rodrigues LR, Galvao RKH, Yoneyama T. System level RUL estimation 
for multiple-component systems. In: Proc 2013 Annu Conf Prognostics Health 
Manage. Soc; 2013. p. 74–82. 
[24] Xing J, Zeng Z, Zio E. A framework for dynamic risk assessment with condition 
monitoring data and inspection data. Reliab Eng Syst Safe 2019;191:106552. 
[25] Wang H, Khan F, Ahmed S, Imtiaz S. Design of Scenario-Based Early Warning 
System for Process Operations. Ind Eng Chem Res 2015;54(33):8255–65. 
[26] Adedigba SA, Oloruntobi O, Khan F, Butt S. Data-driven dynamic risk analysis of 
offshore drilling operations. J Petrol Sci Eng 2018;165:444–52. 
[27] Mamudu A, Khan F, Zendehboudi S, Adedigba S. Dynamic risk assessment of 
reservoir production using data-driven probabilistic approach. J Petrol Sci Eng 
2020;184:106486. 
[28] Zhang X, Zhang L, Hu J. Real-time risk assessment of a fracturing manifold system 
used for shale-gas well hydraulic fracturing activity based on a hybrid Bayesian 
network. J Nat Gas Sci Eng 2019;62:79–91. 
[29] Zywiec WJ, Mazzuchi TA, Sarkani S. Analysis of process criticality accident risk 
using a metamodel-driven Bayesian network. Reliab Eng Syst Safe 2021;207: 
107322. 
[30] Yu Q, Teixeira AP, Liu K, Rong H, Soares CG. An integrated dynamic ship risk 
model based on Bayesian networks and evidential reasoning. Reliab Eng Syst Safe 
2021;216:107993. 
[31] Islam R, Khan F, Venkatesan R. Real time risk analysis of kick detection: Testing 
and validation. Reliab Eng Syst Safe 2017;161:25–37. 
[32] He R, Li X, Chen G, Chen G, Liu Y. Generative adversarial network-based semi- 
supervised learning for real-time risk warning of process industries. Expert Syst 
Appl 2020;150:113244. 
[33] Wang H, Khan F, Ahmed S, Imtiaz S. Dynamic quantitative operational risk 
assessment of chemical processes. Chem Eng Sci 2016;142:62–78. 
[34] Hashemi SJ, Ahmed S, Khan F. Loss functions and their applications in process 
safety assessment. Process Saf Prog 2014;33(3):285–91. 
[35] Sule I, Imtiaz S, Khan F, Butt S. Risk analysis of well blowout scenarios during 
managed pressure drilling operation. J Petrol Sci Eng 2019;182:106296. 
[36] Wang Q, Zhang L, Hu J. Real-time risk assessment of casing-failure incidents in a 
whole fracturing process. Process Saf Environ 2018;120:206–14. 
[37] Dimaio F, Scapinello O, Zio E, Ciarapica C, Cincotta S, Crivellari A, Decarli L, 
Larosa L. Accounting for safety barriers degradation in the risk assessment of oil 
and gas systems by multistate Bayesian networks. Reliab Eng Syst Safe 2021;216: 
107943. 
[38] Wu S, Zhang L, Fan J, Zheng W, Zhou Y. Real-time risk analysis method for 
diagnosis and warning of offshore downhole drilling incident. J Loss Prevent Proc 
2019;62:103933.[39] Parhizkar T, Hogenboom S, Vinnem JE, Utne IB. Data driven approach to risk 
management and decision support for dynamic positioning systems. Reliab Eng 
Syst Safe 2020;201:106964. 
[40] Ahn J, Chang D. Fuzzy-based HAZOP study for process industry. J Hazard Mater 
2016;317:303–11. 
[41] Zadakbar O, Imtiaz S, Khan F. Dynamic risk assessment and fault detection using a 
multivariate technique. Process Saf Prog 2013;32(4):365–75. 
[42] Yu J. A particle filter driven dynamic Gaussian mixture model approach for 
complex process monitoring and fault diagnosis. J Process Contr 2012;22:778–88. 
[43] Figueiredo MAF, Jain AK. Unsupervised learning of finite mixture models. IEEE 
Trans Pattern Anal Mach Intell 2002;24:381–96. 
[44] Arulampalam MS, Maskell S, Gordon N, Clapp T. A tutorial on particle filters for 
online nonlinear/non-Gaussian Bayesian tracking. IEEE Trans Signal Proces 2002; 
50(2):175–88. 
[45] Kotecha JH, Djuric PM. Gaussian sum particle filtering. IEEE Trans Signal Proces 
2003;51(10):2602–12. 
[46] Yu J, Qin SJ. Multimode process monitoring with Bayesian inference-based finite 
Gaussian Mixture models. AIChE J 2008;54(7):1811–29. 
[47] Li C, Mahadevan S, Ling Y, Choze S, Wang L. Dynamic Bayesian network for 
aircraft wing health monitoring digital twin. AIAA J 2017;55(3). 
[48] Ghosh A, Ahmed S, Khan F. Modeling and testing of temporal dependency in the 
failure of a process system. Ind Eng Chem Res 2019;58:8162–71. 
[49] Sule I, Khan F, Butt S, Yang M. Kick control reliability analysis of managed pressure 
drilling operation. J Loss Prevent Proc 2018;52:7–20. 
[50] Abimbola M, Khan F, Khakzad N, Butt S. Safety and risk analysis of managed 
pressure drilling operation using Bayesian network. Safety Sci 2015;76:133–44. 
[51] Yan L, Wu H, Yan Y. Application of fine managed pressure drilling technique in 
complex wells with both blowout and lost circulation risks. Natural Gas Indus B 
2015;2:192–7. 
[52] “Case western reserve university bearing data center website”, https://csegroups. 
case.edu/bearingdatacenter. 
[53] Smith WA, Randall RB. Rolling element bearing diagnostics using the Case Western 
Reserve University data: A bench mark study. Mech Syst Signal Pr 2015;64-65: 
100–31. 
[54] Goyal D, Vanraj Pabla BS, Dhami SS. Condition monitoring parameters for fault 
diagnosis of fixed axis gearbox: A review. Arch Computat Methods Eng 2017;24: 
543–56. 
[55] Karanki D R, Kushwaha H S, Verma A K, Ajit S. Uncertainty analysis based on 
probability bounds (P-box) approach in probabilistic safety assessment. Risk Anal 
2009;29(5):662–75. 
[56] Jones B, Jenkinson I, Yang Z, Wang J. The use of Bayesian network modelling for 
maintenance planning in a manufacturing industry. Reliab Eng Syst Safe 2010;95: 
267–77. 
R. He et al. 
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0011
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0011
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0011
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0012
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0012
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0013
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0013
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0014
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0014
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0015
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0015
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0016
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0016
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0016
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0017
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0017
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0018
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0018
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0019
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0019
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0019
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0020
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0020
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0020
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0021
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0021
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0021
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0022
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0022
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0023
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0023
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0023
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0024
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0024
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0025
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0025
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0026
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0026
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0027
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0027
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0027
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0028
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0028
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0028
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0029
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0029
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0029
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0030
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0030
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0030
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0031
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0031
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0032
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0032
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0032
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0033
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0033
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0034
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0034
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0035
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0035
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0036
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0036
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0037
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0037
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0037
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0037
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0038
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0038
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0038
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0039
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0039
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0039
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0040
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0040
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0041
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0041
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0042
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0042
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0043
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0043
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0044
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0044
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0044
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0045
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0045http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0046
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0046
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0047
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0047
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0048
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0048
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0049
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0049
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0050
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0050
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0051
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0051
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0051
https://csegroups.case.edu/bearingdatacenter
https://csegroups.case.edu/bearingdatacenter
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0053
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0053
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0053
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0054
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0054
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0054
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0055
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0055
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0055
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0056
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0056
http://refhub.elsevier.com/S0951-8320(22)00323-4/sbref0056
	A real-time probabilistic risk assessment method for the petrochemical industry based on data monitoring
	1 Introduction
	2 Preliminaries and problem formulation
	3 Proposed methodology
	3.1 Probability calculation for online signal exceeding alarm threshold
	3.1.1 Probability estimations for environmental nodes
	3.1.2 Probability estimations for process nodes
	3.2 Integrated real-time risk analysis
	4 Implementation and model discussion
	4.1 Case study 1: RT 580 experimental setup
	4.1.1 System description
	4.1.2 PRA results and discussions
	4.2 Case study 2: managed pressure drilling operations
	4.2.1 System description
	4.2.2 Real-time PRA results
	4.2.3 The comparison of GMM and GMM-PF
	4.2.4 Model parameter analysis
	4.2.5 Model uncertainty and validation
	Conclusions
	CRediT authorship contribution statement
	Declaration of Competing Interest
	Acknowledgments
	References