CN104778370A - 基于蒙特卡洛仿真求解动态故障树模型的风险分析方法 - Google Patents

Abstract

本发明公开了一种基于蒙特卡洛仿真求解的动态故障树模型的定量分析方法，所述方法首先对***建立动态故障树模型，而后利用其组件故障概率数据进行故障概率参数动态估计，作为动态故障树运算的输入值，采用蒙特卡洛方法，对所述动态故障树模型进行仿真，通过所述仿真求解动态故障树模型，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。本发明的风险分析方法不仅可以对状态数量巨大的复杂***进行有效的风险估计，同时也可以得到***组件的重要度分析数据，即同时也可以对组件的故障率未知及其故障率的变化情况进行有效的风险分析。

Description

基于蒙特卡洛仿真求解动态故障树模型的风险分析方法

技术领域

本发明属于***安全及风险分析技术领域，具体涉及一种基于蒙特卡洛仿真求解动态故障树的风险分析方法。

背景技术

随着计算机技术的不断普及，各类控制***在各行各业发挥着越来越重要的作用，如飞行控制***、列车控制***、生产控制***等。由于控制***的重要性，其***的安全可靠性就显得尤为重要，而准确的分析***面临的风险，对***的开发和安全运行具有重要意义。

例如，得益于计算机技术、电子技术和控制技术的飞速发展，列车控制***获得了很好的发展机遇，目前已经在高速铁路上广泛使用CTCS-3基于通信的列车控制***，其主要包括三个模块，即地面子***(轨道电路、列控中心、车站计算机联锁、列控中心等)、车载子***(安全计算机、测速定位、无线通信模块、人—机接口等)和信息传输子***(如GSM—R等)。在铁路不断向着高速、高密度的方向发展的现实前提下，对控制***这样由一系列子***模块组成来控制列车运行速度，为了保障列车运行安全和提高运输能力，列车控制***需要满足的安全条件愈发严酷，控制要求也越来越高，***的结构也日趋复杂。

随着***包含的组件也越来越多的情况下，虽然列车控制***在保障列车安全可靠高速运行方面有着优异的表现，但是存在的风险也不容忽视。由于***中组件数量巨大且***使用数量亦较大，由组件的随机故障导致的***风险也不断加大。列控***一旦发生故障，将会带来巨大的经济损失，且有可能危机乘客的生命安全，造成严重的社会影响。

目前，对控制***的风险分析方法有几十种之多，近年来，风险分析方法向正向定性与定量结合、注重动态***描述能力、计算机仿真等方向发展。在现有技术中，对列车控制***的风险分析，在利用了概率方法和衡量指标之后得到广泛的发展，方法和指标的选择需要根据具体的问题及假设来选定，而风险分析所使用的模型有效性，直接关系到风险分析的有效性。在选用合适的风险分析模型时，也要考虑到实际的故障分布无法与分析中所用的分布模式完全匹配的问题，从而避免在简化时引入关键性错误。因此，目前常采用动态故障树模型对列车控制***进行风险分析。

在对于动态故障树模型的求解上，通常采用将动态故障树模型转换为马尔科夫模型从而使用解析法进行求解。但当动态故障树规模增大时，马尔科夫模型中的状态量增加速度太快，存在状态***问题。虽然以往的研究者提出了各种改进方法，但对于存在大量状态的复杂***仍然应用有限。另外，将动态故障树模型转化为马尔科夫模型，使用解析法来求解时，对于组件的故障率均只能假设为常数，而实际的控制***运行过程中组件的故障率往往随着时间的变化有所变化，这就使得所求解的组件的故障率偏离真实值，无法真实的反映组件故障率的变化情况，导致无法进行有效的风险估计。

发明内容

本发明的目的是提供一种基于蒙特卡洛仿真求解动态故障树模型的风险分析方法，对状态数量巨大的复杂***进行有效的风险估计，同时得到组件的重要度分析数据，对组件的故障率及其故障率的变化情况进行有效的风险分析。

根据本发明的一个方面，提供了一种基于蒙特卡洛仿真求解动态故障树模型的风险分析方法，所述方法包括：

对需要进行风险分析的***建立动态故障树模型；

估计所述***的组件故障概率参数；

将所述组件故障概率参数作为动态故障树模型运算的输入值，采用蒙特卡洛仿真方法对所述动态故障树模型进行求解，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。

上述方案中，所述建立动态故障树模型，进一步包括：根据所述***的结构，以***故障为顶事件，以所述***的组件故障为基事件，建立所述***的动态故障树模型。

上述方案中，所述估计所述***的组件故障概率参数，进一步包括：根据组件的故障时间数据，采用预设的参数估计，得到组件故障概率参数。

上述方案中，所述预设的参数估计，包括最小二乘估计、最大似然估计、贝叶斯估计中的一种或多种。

上述方案中，所述预设的参数估计为贝叶斯最大似然估计。

上述方案中，所述贝叶斯最大似然估计的参数估计过程中，通过马尔科夫-蒙特卡洛方法建立马尔科夫链，使马尔科夫链的稳定分布和贝叶斯最大似然估计的后验分布一致，当马尔科夫链收敛时，将马尔科夫链的模拟值作为从后验分布中抽取的样本。

上述方案中，所述采用蒙特卡洛仿真方法对所述动态故障树模型进行求解，进一步包括：

选取仿真时间步长；

依据概率分布或概率密度判断组件状态；

根据所述组件状态对所述动态故障树进行逻辑门运算。

上述方案中，所述时间步长为固定步长或可变步长。

上述方案中，所述依据概率分布判断组件状态，进一步包括：

将仿真时间与所述组件的故障时间进行比较，判断当前所述时间步长内的组件状态。

上述方案中，所述依据概率密度判断组件状态，进一步包括：

设定所述组件的故障概率为f；

根据所述组件的故障概率密度函数，在预设时间区间上抽取服从均匀分布的随机数r，若r≤f，则判断所述组件故障；若r＞f，则判断所述组件非故障。

本发明实施例的基于蒙特卡洛仿真求解的动态故障树模型的风险分析方法，首先对需要进行风险分析的***建立动态故障树模型，而后估计所述***的组件故障概率数据，并将所述组件故障概率数据作为动态故障树运算的输入值，采用蒙特卡洛方法，对所述动态故障树模型进行仿真，通过所述仿真求解动态故障树模型，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。本发明的风险分析方法不仅可以对状态数量巨大的复杂***进行有效的风险估计，同时也可以得到***组件的重要度分析数据，即同时也可以对组件的故障率及其故障率的变化情况进行有效的风险分析。

附图说明

为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明第一实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图；

图2为图1所示风险分析方法的流程框图；

图3为本发明第二实施例的基于蒙特卡洛方仿真求解动态故障树模型的风险分析方法流程示意图；

图4a为图3所示风险分析方法的参数估计过程的α值随迭代次数的变化图；

图4b为图3所示风险分析方法的参数估计过程的实现中的α取值直方图；

图4c为图3所示风险分析方法的参数估计过程的β值随迭代次数的变化图；

图4d为图3所示风险分析方法的参数估计过程的实现中的β取值直方图；

图5为图3所示的风险分析方法的参数估计中最大似然估计和贝叶斯最大似然估计比较图；

图6为图3所示风险分析方法的仿真流程示意图；

图7为本发明第三实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图；

图8为图7所示风险分析方法的仿真流程示意图；

图9a为图7所示风险分析方法的逻辑门计算中的优先与门时序状态图；

图9b为图7所示风险分析方法的逻辑门计算中的热备门时序状态图；

图9c为图7所示风险分析方法的逻辑门计算中的顺序门时序状态图；

图9d为图7所示风险分析方法的逻辑门计算中的功能相关门时序状态图；

图10为本发明第四实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图；

图11为图10所示风险分析方法的仿真流程示意图；

图12a为本发明第三实施例采用概率分布抽取法时的仿真流程示意图；

图12b为本发明第四实施例采用概率密度抽取法时的仿真流程示意图。

具体实施方式

本技术领域技术人员可以理解，除非特意声明，这里使用的单数形式“一”、“一个”、“所述”和“该”也可包括复数形式。应该进一步理解的是，本发明的说明书中使用的措辞“包括”是指存在所述特征、整数、步骤、操作、元件和/或组件，但是并不排除存在或添加一个或多个其他特征、整数、步骤、操作、元件、组件和/或它们的组。应该理解，当我们称元件被“连接”或“耦接”到另一元件时，它可以直接连接或耦接到其他元件，或者也可以存在中间元件。此外，这里使用的“连接”或“耦接”可以包括无线连接或耦接。这里使用的措辞“和/或”包括一个或更多个相关联的列出项的任一单元和全部组合。

本技术领域技术人员可以理解，除非另外定义，这里使用的所有术语(包括技术术语和科学术语)具有与本发明所属领域中的普通技术人员的一般理解相同的意义。还应该理解的是，诸如通用字典中定义的那些术语应该被理解为具有与现有技术的上下文中的意义一致的意义，并且除非像这里一样定义，不会用理想化或过于正式的含义来解释。

为便于对本发明实施例的理解，下面详细描述本发明的实施方式，通过参考附图描述的实施方式是示例性的，仅用于解释本发明，而不能解释为对本发明的限制。

本发明针对控制***的安全问题对风险估计方法进行优化，尤其针对列车控制***，采用蒙特卡洛方法对控制***的求解模型-动态故障树模型进行求解，从而得到相应***的故障概率分布分析数据，同时得到***的组件风险分析数据，一般指的是***组件的重要度分析数据。下面结合具体实施例及附图，对本发明作进一步的详细说明。

图1为本发明第一实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图。

如图1所示，本实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法，包括如下步骤：

步骤S11，对需要进行风险分析的***建立动态故障树模型。

优选的，本步骤中，构建动态故障树模型，根据所述***的结构，以***故障为顶事件，以所述***的组件故障为基事件，建立所述***的动态故障树模型。

步骤S12，估计所述***的组件故障概率参数。

优选的，本步骤中，根据组件的故障时间数据，采用预设的参数估计，得到组件故障概率数据。所述预设的参数估计，包括最小二乘估计、最大似然估计、贝叶斯估计中的一种或多种。进一步的，优选为贝叶斯最大似然估计。这里的贝叶斯最大似然估计，是最大似然估计与贝叶斯估计相结合得到的一种新的参数估计方法，兼具贝叶斯估计和最大似然估计的优点。

优选的，所述贝叶斯最大似然估计的参数估计过程中，通过马尔科夫-蒙特卡洛方法建立马尔科夫链，使马尔科夫链的稳定分布和贝叶斯最大似然估计的后验分布一致，当马尔科夫链收敛时，将马尔科夫链的模拟值作为从后验分布中抽取的样本。

步骤S13，将所述组件故障概率参数作为动态故障树模型运算的输入值，采用蒙特卡洛仿真方法，对所述动态故障树模型进行求解。

本步骤中，通过蒙特卡洛仿真的方法，对所建立的动态故障树进行求解，不仅可以避免出现状态***的情况，同时也可以对组件进行故障概率计算，即组件的重要度分析参数。

对于动态故障树的蒙特卡洛仿真，按照仿真的时间步长分，有固定时间步长和可变时间步长两种方式来表示模拟中的时间进程。按照仿真中对于组件状态即基事件状态的判断方法分，可分为按照组件故障概率分布抽取和按照组件的概率密度抽取两种方法。

步骤S14，通过所述仿真求解动态故障树模型，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。

当对列车控制(简称“列控”)***采用本实施例进行风险估计时，通过上述过程，可以得到列控***组件重要度分析数据，即列控***组件的灵敏度分析。

在列车控制***风险分析中，特定组件在***的架构、位置可能造成其影响***安全可靠运行的程度，所以分析组件的重要度在***设计、诊断及最优化分析时有相当大的作用，在***检查、维护及故障检测执行的先后顺序可根据组件重要度分析出的各组件灵敏度作为依据，或是在***改进时改进重要度较大的组件。

灵敏度分析大致可有三种方法：(1)数学模型计算法、(2)统计灵敏度分析法、(3)图形灵敏度分析法。优选的，本实施例采用精确度相对较高的数学模型计算法。

数学模型计算法评估输入参数在一定范围内变动时，输出的灵敏程度，此方法通常计算输入参数的某些可能值范围，其输出的响应。

本实施例对***风险的分析通过蒙特卡洛仿真来实现，因此在对组件的灵敏度分析时亦采用蒙特卡洛仿真方法，定义***发生故障的时间对组件发生故障的时间灵敏度为：

$S_f=\frac{\mathrm{\ΔE}\left(t_{\mathrm{sys}}\right)}{\mathrm{\ΔE}\left(t_c\right)}$

其中E(t_sys)表示***发生故障的时间期望，E(t_c)表示组件发生故障的时间期望。

上式实质上是***发生故障的时间期望对于组件发生故障的时间期望的偏微分，可以反映出组件故障对***发生故障的时间期望的影响程度。通过该指标可以得到组件对***故障贡献大小的顺序，以此作为寻找提高***可靠性方法的参考。

从以上分析可以看出，本实施例的基于蒙特卡洛方法求解动态故障树模型的风险估计方法，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。本发明的风险分析方法不仅可以对状态数量巨大的复杂***进行有效的风险估计，同时也可以得到***组件的重要度分析数据，即同时也可以对组件的故障率及其故障率的变化情况进行有效的风险分析。

图2为图1所示风险分析方法的流程框图。

如图2所示，对图1所示风险分析方法做了详细说明。这里，进一步以列车控制***为例，对图2进行说明。本实施例的基于动态故障树和蒙特卡洛仿真的列车控制***风险分析是对列车控制***进行动态故障树仿真以得到***故障概率分布及组件重要度的方法。具体做法是：在了解列车控制***的***结构的基础上，对***建立以***故障为顶事件，以***组件故障为基事件的动态故障树。然后利用***运行中组件故障数据的统计，通过参数估计方法得到组件的故障概率分布。接着对动态故障树进行蒙特卡洛仿真，仿真结果经过统计分析得出***的故障概率分布。最后使用改变组件的概率分布函数的动态故障树蒙特卡洛仿真得到各个组件的重要度。

图3为本发明第二实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图。

如图3所示，本实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法，包括如下步骤：

步骤S21，对需要进行风险分析的***建立动态故障树模型。

本步骤与第一实施例构建动态故障树模型的过程相同，在此不再赘述。

步骤S22，通过贝叶斯最大似然估计对组件故障概率参数进行参数估计。

本实施例利用马尔科夫蒙特卡洛(MCMC)方法，结合最大似然估计的两参数威布尔分布参数估计方法。完全样本下的参数贝叶斯最大似然估计如下：

假定同一组件的故障分布相互独立且均服从同一威布尔概率分布，其生存函数为：

$\stackrel{\‾}{F}\left(x\right)=P(X>x)=\exp [-{\left(\frac{x}{\mathrm{\β}}\right)}^{\mathrm{\α}}]---\left(1\right)$

式(1)中x,α,β>0。x的概率密度函数为：

$f\left(x\right)=\mathrm{\α}{\mathrm{\β}}^{-\mathrm{\α}}{x}^{\mathrm{\α}-1}\exp [-{\left(\frac{x}{\mathrm{\β}}\right)}^{\mathrm{\α}}]---\left(2\right)$

式(2)中，f(x)表示组件运行到时间x时的故障率，r(x)表示组件在运行到时间x时依然工作的概率，定义为在时间x上组件的故障密度的条件概率，亦即到时间x时组件仍未故障。因而：

$r\left(x\right)=\frac{f\left(x\right)}{\stackrel{\‾}{F}\left(x\right)}=\mathrm{\α}{\mathrm{\β}}^{-\mathrm{\α}}{x}^{\mathrm{\α}-1}---\left(3\right)$

式(3)中，利用组件的故障时间样本数据(x₁,x₂,……,x_n)来估计分布模型的参数。令θ时表示分布模型的参数,Π(θ)表示先验概率模型的密度函数。根据贝叶斯公式，样本数据和模型参数的联合分布密度函数为：

π(x,θ)＝π(x|θ)π(θ)＝π(θ|x)π(x)

模型参数的后验模型可表示为π(θ|x),则

π(θ|x)＝π(x|θ)π(θ)/π(x)

又因为π(x)和θ相互独立，所以可以认为

π(θ|x)∝π(x|θ)π(θ)

其中π(x|θ)为以θ为条件的似然函数。

但是由于贝叶斯方法要求基于参数的后验分布来进行统计推断，而直接处理后验分布非常困难。所以使用了马尔科夫蒙特卡洛(MCMC)方法——用此方法建立一个马尔科夫链，使它的稳定分布和后验分布一致，当马尔科夫链收敛时的模拟值可以看作从后验分布中抽取的样本。

在以下的说明中，仍以列车控制***为例进行说明。以下为根据某列控组件发生故障的时间样本，表1为列控组件发生故障的时间样本，所述时间样本的单位为小时(h)。如表1所示，对该类组件的故障概率分布参数估计为例介绍参数估计方法的实现：

表1

有些组件的故障率和时间无关(α＝1)，但是对多数组件而言，由于日久损耗或者是其他原因导致，故障率是随着时间增加的(α>1)。

若一套组件的故障时间数据未知，通常会认为故障率是随着时间的增加而增加的(α>1)，极限情形是呈线性关系(α＝2)。此外，α最可能的值近似于1.5，先验概率在1.5两侧以同样的速度递减。因而，α的先验边缘概率估计可选为：

$g\left(\mathrm{\&alpha;}\right)=\left\{\begin{array}{c}4(\mathrm{\&alpha;}-1),1<\mathrm{\&alpha;}<1.5\\ 4(2-\mathrm{\&alpha;}),1.5<\mathrm{\&alpha;}\&le;2\end{array}\right.---\left(4\right)$

这是(1,2)上的对称三角分度概率密度。从这个概率密度上采样，可使R₁,R₂～U(0,1)，取

$\mathrm{\&alpha;}=1+\frac{1}{2}(R_1+R_2)---\left(5\right)$

即为参数α的抽样公式。根据组件发生故障的时间统计，可认为组件发生故障的事件期望在2000h到3000h之间某处。因此，给出

$E\left(X\right|\mathrm{\&alpha;},\mathrm{\&beta;})=\underset{0}{\overset{\&infin;}{\&Integral;}}F\left(x\right)\mathrm{dx}=\mathrm{\&beta;\&Gamma;}(\frac{1}{\mathrm{\&alpha;}}+1)---\left(6\right)$

由上式可以给出β关于α的条件先验概率，即

U(2000/Γ(1/α+1),3000/Γ(1/α+1))分布密度。一旦α已经抽取，β的抽样即可根据下式抽出。

$\mathrm{\&beta;}=\frac{1000(2+R_3)}{\mathrm{\&Gamma;}(1/\mathrm{\&alpha;}+1)}---\left(7\right)$

其中，R₃～U(0,1)。

对参数的贝叶斯估计构造MCMC方法:

α和β的联合先验概率为：

$g(\mathrm{\&alpha;},\mathrm{\&beta;})=\left\{\begin{array}{c}\frac{4}{1000}(\mathrm{\&alpha;}-1)\mathrm{\&Gamma;}(\frac{1}{\mathrm{\&alpha;}}+1),1<\mathrm{\&alpha;}<1.5\\ \frac{4}{1000}(2-\mathrm{\&alpha;})\mathrm{\&Gamma;}(\frac{1}{\mathrm{\&alpha;}}+1),1.5<\mathrm{\&alpha;}\&le;2\end{array}\right.---\left(8\right)$

其中，2000/Γ(1/α+1)＜β＜3000/Γ(1/α+1)。

根据贝叶斯估计知，α和β的后验概率为π(α,β)∝L(α,β)g(α,β)，其中L(α,β)是最大似然函数。

为了得到参数α和β的估计值，用MCMC方法根据π(α,β)对α和β采样，其建议分布取α和β的联合分布g(α,β)。当马尔科夫链接***衡条件时便可被取样。

接受概率为：

$\min =[1,\frac{\mathrm{\&pi;}({\mathrm{\&alpha;}}^{\&prime;},{\mathrm{\&beta;}}^{\&prime;})g(\mathrm{\&alpha;},\mathrm{\&beta;})}{\mathrm{\&pi;}(\mathrm{\&alpha;},\mathrm{\&beta;})g({\mathrm{\&alpha;}}^{\&prime;},{\mathrm{\&beta;}}^{\&prime;})}]=\min [1,\frac{L({\mathrm{\&alpha;}}^{\&prime;},{\mathrm{\&beta;}}^{\&prime;})}{L(\mathrm{\&alpha;},\mathrm{\&beta;})}]---\left(9\right)$

其中，

$\begin{array}{c}L(\mathrm{\&alpha;},\mathrm{\&beta;})=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}\mathrm{\&alpha;}{x_i}^{\mathrm{\&alpha;}-1}{\mathrm{\&beta;}}^{-\mathrm{\&alpha;}}\exp [-{\left(\frac{x_i}{\mathrm{\&beta;}}\right)}^{\mathrm{\&alpha;}}]\\ ={\mathrm{\&alpha;}}^n{\mathrm{\&beta;}}^{-\mathrm{n\&alpha;}}\exp [\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}-{\left(\frac{x_i}{\mathrm{\&beta;}}\right)}^{\mathrm{\&alpha;}}]{(x_1...x_n)}^{\mathrm{\&alpha;}-1}\end{array}---\left(10\right)$

在实现中步骤为：

步骤A，给定参数α,β的初始值，根据α,β的先验概率抽取α′,β′；

步骤B，计算接受概率；

步骤C，以接受概率接受抽取的参数值，反之，则拒绝备选参数值；

步骤D，重复步骤A～C，直到收敛。

图4a为图3所示风险分析方法的参数估计过程的α值随迭代次数的变化图；图4b为图3所示风险分析方法的参数估计过程的实现中的α取值直方图；图4c为图3所示风险分析方法的参数估计过程的β值随迭代次数的变化图；图4d为图3所示风险分析方法的参数估计过程的实现中的β取值直方图。如图4a至图4d所示，此组件的故障概率分布贝叶斯参数估计迭代实现：取迭代次数为5000，实现以上算法，其中α值和β值随着迭代次数的变化图分别如图4a和图4b所示。α值和β值的估计值分别取迭代实现中α和β所有取值的平均值，得α为1.14，β为2495.1。图4c和图4d分别为迭代实现中α取值的直方图和β取值的直方图。

图5为图3所示的风险分析方法的参数估计中最大似然估计和贝叶斯最大似然估计比较图。如图5所示，其中，柱状图代表最大似然估计法所得的组件故障概率密度函数变化情况，曲线代表贝叶斯最大似然估计所得的组件故障概率密度函数曲线。使用最大似然估计法对此样本进行的参数估计计算得α值为1，β值2201。由先验知识可知该组件的故障率随着时间的变化而变化，而由最大似然估计法得到的α值为1表示组件的故障率是一常数，不随时间变化。故而相比之下贝叶斯最大似然估计法所得参数值较为准确。另从组件故障时间样本的直方图和两种估计法的故障概率密度曲线直接观察亦可知贝叶斯最大似然估计更为精确。

这里需要说明的是，在参数估计的具体算法实现中，为了计算便利，将接受概率做了变换，使其为lnL(α′,β′)-lnL(α,β)。

优选的，本实施例还针对截尾数据样本进行了参数贝叶斯估计，具体估计过程如下：

对于寿命周期较长的安全相关***中组件，由于设备的定期维修或者更换等原因，可收集到的故障或失效数据在整体上往往更符合随机截尾数据类型的范围。所以在完全样本的参数估计之外讨论截尾数据样本下的参数估计。

假设随机截尾的数据样本所表示的组件故障概率分布为F(x,θ)，分布密度为f(x,θ)，其中θ为未知参数向量。可将随机截尾的数据样本整理为(t₁,δ₁)，(t₂,δ₂)，……，(t_n,δ_n)。其中δ_i＝1表示寿终数据，δ_i＝0表示截尾数据。t_i为x变量的取值，δ_i为θ变量的取值，i为1到n的自然数。

这时的似然函数是:

$L\left(\mathrm{\&theta;}\right)=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}{\left(f(t_i,\mathrm{\&theta;})\right)}^{{\mathrm{\&delta;}}_i}{[1-F(t_i,\mathrm{\&theta;})]}^{1-{\mathrm{\&delta;}}_i}---\left(11\right)$

在对服从威布尔分布的组件故障率随机截尾样本的参数估计上，依然适用上述完全样本参数估计中基于MCMC的贝叶斯最大似然估计方法，仅改变其最大似然估计函数为：

$L(\mathrm{\&alpha;},\mathrm{\&beta;})=\underset{i=1}{\overset{n}{\mathrm{\&Pi;}}}{\left\{\mathrm{\&alpha;}{x_i}^{\mathrm{\&alpha;}-1}{\mathrm{\&beta;}}^{-\mathrm{\&alpha;}}\exp \right[-{\left(\frac{x_i}{\mathrm{\&beta;}}\right)}^{\mathrm{\&alpha;}}\left]\right\}}^{{\mathrm{\&delta;}}_i}\&CenterDot;{\left\{\exp \right[-{\left(\frac{x_i}{\mathrm{\&beta;}}\right)}^{\mathrm{\&alpha;}}\left]\right\}}^{1-{\mathrm{\&delta;}}_i}---\left(12\right)$

式(12)中，α是形状参数，β是比例参数。

步骤S23，选取固定时间步长为仿真推进时间。

步骤S24，依据概率分布判断组件状态。

本步骤中，依据概率分布判断组件状态，进一步包括：将仿真时间与所述组件的故障时间进行比较，判断当前所述时间步长内的组件状态。

步骤S25，根据所述组件状态对所述动态故障树进行逻辑门运算。

步骤S26，通过上述仿真求解动态故障树模型，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。

图6为图3所示风险分析方法的仿真流程示意图。

如图6所示，本实施例的风险分析方法，采用固定时间步长，这里对仿真具体过程的说明，仍以列车控制***为例。

仿真固定时间步长法模拟时钟是按照一个固定的仿真时间片作为时间增量来推进仿真，在每个时间片推进点上对各个事件的状态进行评估，更新动态故障树的仿真状态。

采用组件的概率分布抽取法是指在按照组件的概率分布函数来抽取组件发生故障的时间，在仿真中以仿真时间是否大于组件发生故障的时间来判断组件的状态，作为动态故障树运算的输入。

固定步长动态故障树蒙特卡洛仿真步骤如下：

设仿真总次数M，仿真时钟时间片T，用m表示仿真次数，t表示仿真时间，t_f1,t_f2,···t_fn表示组件发生故障的时间，f_s＝1表示***故障。

仿真开始，置m＝m+1。判断m是否不小于M，是则结束仿真。

按照各个组件的故障概率分布来抽取组件发生故障的时间亦即基事件发生的时间t_f1,t_f2,···t_fn。

随着时间按照仿真时间片T为单位的推进，根据各个基事件状态按照逻辑门运算来更新***动态故障树状态，确定***状态。

统计***在M次仿真中发生故障的时间。使用先验分布来估计***的故障概率分布。

图7为本发明第三实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图。

如图7所示，本实施例的风险分析方法，包括如下步骤：

步骤S31，对需要进行风险分析的***建立动态故障树模型。

本步骤与第一实施例构建动态故障树模型的过程相同，在此不再赘述。

步骤S32，通过贝叶斯最大似然估计对组件故障概率数据进行参数估计。

步骤S33，选取可变时间步长为仿真推进时间。

本步骤与第二实施例的基本相同，所不同的是，本实施例在仿真过程中所选取的仿真推进时间步长是可变时间步长。

图8为图7所示风险分析方法的仿真流程示意图。

如图8所示，本实施例的风险分析方法，采用可变时间步长，这里对仿真具体过程的说明，仍以列车控制***为例。

可变步长法仿真中，仿真的时间增量由动态故障树中事件出现的时间间隔确定，即把事件出现的时刻作为仿真时钟的推进点，并在该时刻后对事件及***状态进行评估和更新故障树状态。

可变步长动态故障树蒙特卡洛仿真的步骤如下：

仿真开始，置m＝m+1。判断m是否不小于M，是则结束仿真。

按照各个组件的故障概率分布来抽取组件发生故障的时间亦即基事件发生的时间t_f1,t_f2,……t_fn。

对t_f1,t_f2,……t_fn进行排序，最小时间为t_min，根据各个基事件状态按照逻辑门运算来更新***动态故障树状态，确定***状态，并按照Δt推进***仿真时间。Δt为可变时间步长，取值为经排序后相邻两个故障时间的差值。

统计***在M次仿真中发生故障的时间。使用先验分布来估计***的故障概率分布。

步骤S34，从所述组件故障概率数据中分布式抽取故障概率计算所述仿真时间步长内的组件故障率，并根据所述组件故障率判断组件状态。

步骤S35，根据所述组件状态对所述动态故障树进行逻辑门运算。

在仿真实现步骤中的动态故障树逻辑门运算是根据各个基事件的状态即组件故障是否发生来判断顶事件***故障是否发生的运算过程。其中关于静态逻辑门的运算简单，只用进行布尔值运算就可以完成，但是动态逻辑门的运算由于与时序逻辑有关，所以需要重点考虑。以下结合附图9a至附图9d，从动态逻辑门的时序逻辑出发对其在态故障树蒙特卡洛仿真中的实现进行描述。

图9a为图7所示风险分析方法的逻辑门计算中的优先与门时序状态图。

如图9a所示，优先与门当其基事件按照预先设定的顺序(通常为从左到右)依次发生，即组件按顺序依次故障之后才会导致输出端事件的发生。如下图所示，在场景1和场景2中，组件A和组件B依次故障，导致了输出端事件的发生。而场景3中组件B先于组件A故障，因此输出端事件不发生。在蒙特卡洛仿真实现中，某时刻t时优先与门的输出既需要考虑该时刻与门基事件的状态，还需要考虑前一仿真时刻t′的与门基事件的状态，仅当组件的故障态依次出现时输出端事件才发生。例如下图假设基事件为A和B，用0和1分别正常状态和故障状态，则仅当AB两事件的状态从10转到11时优先与门输出事件发生。

图9b为图7所示风险分析方法的逻辑门计算中的热备门时序状态图。

如图9b所示，热备门中分别存在主件即工作件和热备件两部分。主件备件同时处于运行状态，但只有主件作为输出，备件在主件失效时启用作为输出。但是备件在主件在工作态时亦存在失效的可能性。于是，热备门中的输出事件发生可以在两种情况下发生：一是主件失效转换至备件工作后，某个时刻在主件并未修复的情况下备件亦发生失效；二是主件失效转换至备件时备件已经失效。两种失效情形如下图所示。在蒙特卡洛仿真实现中热备门在仿真时刻t的输出事件状态即可根据上述两种情况判断。

图9c为图7所示风险分析方法的逻辑门计算中的顺序门时序状态图。

如图9c所示，顺序门的失效方式与优先与门类似，但是顺序门中事件的发生必须按照特定的方式实现，即首组件的失效导致次组件顶替进入工作状态，接着依次后移，直到顺序门所有组件都失效，则输出事件发生。没有组件能先于首组件失效。其失效情形如下图所示。若在组件可修复的情况下，若在所有组件都失效导致顺序门输出事件发生的情况下某个组件被修复，则***转换到工作态。在蒙特卡洛仿真中，在首组件失效时才进行次组件的时间计算，依次直到所有组件都失效时，才改变顺序门的输出。

图9d为图7所示风险分析方法的逻辑门计算中的功能相关门时序状态图。

如图9d所示，功能相关门的表示触发事件和其他依赖事件的逻辑关系，当触发事件发生时，功能相关门关联的依赖事件都会发生。在触发事件发生时，依赖事件都会位于实际上的失效状态。

步骤S36，通过上述仿真求解动态故障树模型，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。

图10为本发明第四实施例的基于蒙特卡洛仿真求解动态故障树模型的风险分析方法流程示意图。

如图10所示，本实施例的风险分析方法包括如下步骤：

步骤S41，对需要进行风险分析的***建立动态故障树模型。

本步骤与第一实施例构建动态故障树模型的过程相同，在此不再赘述。

步骤S42，通过贝叶斯最大似然估计对组件故障概率数据进行参数估计。

步骤S43，选取可变时间步长为仿真推进时间。

步骤S44，依据概率密度判断组件状态。

本步骤中，所述依据概率密度判断组件状态，进一步包括：设定所述组件的故障概率为f；根据所述组件的故障概率密度函数，在预设时间区间上抽取服从均匀分布的随机数r，若r≤f，则判断所述组件故障；若r＞f，则判断所述组件非故障。

按照组件的概率密度抽取法是指在仿真中，某个仿真时间t上，按照组件的概率密度得到该时刻组件的故障概率，然后在(0,1)区间抽取服从平均分布的随机数来判断该时刻组件的状态，作为动态故障树运算的输入。

图11为图10所示风险分析方法的仿真流程示意图；图12a为本发明第三实施例采用概率分布抽取法时的仿真流程示意图；图12b为本发明第四实施例采用概率密度抽取法时的仿真流程示意图。

如图11所示，按照组件的故障概率密度抽取法的不同在于不抽取组件发生故障的时间，而是在仿真中按照时刻t组件的故障概率来确定组件的状态。假设根据组件的故障概率密度函数，时刻t该组件的故障率为f，则抽取在(0,1)上抽取服从均匀分布的随机数r，若r≤f，则该组件故障，否则该组件非故障。两种方法在实现步骤中的区别如图12a和图12b所示。从图12a可以看出，采用概率分布抽取法时，首先按照概率分布确定基事件的发生时间，在仿真进行中当仿真时间到达基事件时间后认为组件状态为故障，否则认为组件正常工作；从图12b可以看出，采用概率密度抽取法时，需要在仿真进行当中(即每个仿真步长内)，计算该时刻组件的故障率并随机判断组件是否故障。

步骤S45，根据所述组件状态对所述动态故障树进行逻辑门运算。

步骤S46，通过上述仿真求解动态故障树模型，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。

本实施例的基于蒙特卡洛方法求解动态故障树模型的风险估计方法，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。本发明的风险分析方法为状态空间巨大的复杂***的定量风险分析提供了一种高效自动求解方法。与其他同类算法相比，有效的扩充了状态空间的分析能力，提高了分析效率。同时，本发明的风险分析方法可以利用组件的可靠性统计数据或在线数据进行自动分析，可以解决组件可靠性参数未知或由于可靠性参数动态变化无法解析表达这一困扰大部分新型设备定量风险分析的问题。

本领域普通技术人员可以理解：附图只是一个实施例的示意图，附图中的模块或流程并不一定是实施本发明所必须的。

通过以上的实施方式的描述可知，本领域的技术人员可以清楚地了解到本发明可借助软件加必需的通用硬件平台的方式来实现。基于这样的理解，本发明的技术方案本质上或者说对现有技术做出贡献的部分可以以软件产品的形式体现出来，该计算机软件产品可以存储在存储介质中，如ROM/RAM、磁碟、光盘等，包括若干指令用以使得一台计算机设备(可以是个人计算机，服务器，或者网络设备等)执行本发明各个实施例或者实施例的某些部分所述的方法。

本说明书中的各个实施例均采用递进的方式描述，各个实施例之间相同相似的部分互相参见即可，每个实施例重点说明的都是与其他实施例的不同之处。尤其，对于装置或***实施例而言，由于其基本相似于方法实施例，所以描述得比较简单，相关之处参见方法实施例的部分说明即可。以上所描述的装置及***实施例仅仅是示意性的，其中所述作为分离部件说明的单元可以是或者也可以不是物理上分开的，作为单元显示的部件可以是或者也可以不是物理单元，即可以位于一个地方，或者也可以分布到多个网络单元上。可以根据实际的需要选择其中的部分或者全部模块来实现本实施例方案的目的。本领域普通技术人员在不付出创造性劳动的情况下，即可以理解并实施。

以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求的保护范围为准。

Claims (10)

1.一种基于蒙特卡洛仿真求解动态故障树模型的风险分析方法，其特征在于，所述方法包括：

对需要进行风险分析的***建立动态故障树模型；

估计所述***的组件故障概率参数；

将所述组件故障概率参数作为动态故障树模型运算的输入值，采用蒙特卡洛仿真方法对所述动态故障树模型进行求解，得到所述***的故障概率分布分析数据和所述***的组件重要度分析数据。

2.根据权利要求1所述的风险分析方法，其特征在于，所述建立动态故障树模型，进一步包括：根据所述***的结构，以***故障为顶事件，以所述***的组件故障为基事件，建立所述***的动态故障树模型。

3.根据权利要求1所述的风险分析方法，其特征在于，所述估计所述***的组件故障概率参数，进一步包括：根据组件的故障时间数据，采用预设的参数估计得到组件故障概率参数。

4.根据权利要求3所述的风险分析方法，其特征在于，所述预设的参数估计包括最小二乘估计、最大似然估计、贝叶斯估计中的一种或多种。

5.根据权利要求3所述的风险分析方法，其特征在于，所述预设的参数估计为贝叶斯最大似然估计。

6.根据权利要求5所述的风险分析方法，其特征在于，所述贝叶斯最大似然估计的参数估计过程中，通过马尔科夫-蒙特卡洛方法建立马尔科夫链，使马尔科夫链的稳定分布和贝叶斯最大似然估计的后验分布一致，当马尔科夫链收敛时，将马尔科夫链的模拟值作为从后验分布中抽取的样本。

7.根据权利要求1所述的风险分析方法，其特征在于，所述采用蒙特卡洛仿真方法对所述动态故障树模型进行求解，进一步包括：

选取仿真时间步长；

依据概率分布或概率密度判断组件状态；

根据所述组件状态对所述动态故障树进行逻辑门运算。

8.根据权利要求7所述的风险分析方法，其特征在于，所述时间步长为固定步长或可变步长。

9.根据权利要求7所述的风险分析方法，其特征在于，所述依据概率分布判断组件状态，进一步包括：

将仿真时间与所述组件的故障时间进行比较，判断当前所述时间步长内的组件状态。

10.根据权利要求7所述的风险分析方法，其特征在于，所述依据概率密度判断组件状态，进一步包括：

设定所述组件的故障概率为f；

根据所述组件的故障概率密度函数，在预设时间区间上抽取服从均匀分布的随机数r，若r≤f，则判断所述组件故障；若r＞f，则判断所述组件非故障。