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
