陈 妍，何 亮，余少青，陈海英，张春明

(环境保护部核与辐射安全中心，北京100082)



应用分层模型进行核电站设备可靠性参数估计

陈 妍，何 亮，余少青，陈海英，张春明

(环境保护部核与辐射安全中心，北京100082)

目前国内对于从多个核电站统计的设备失效数据进行各特定核电站设备可靠性参数估计的方法研究尚少。本文研究了用于可靠性参数估计的分层模型以及实现分层模型的两种方法：带Kass-Steffey修正的参数经验贝叶斯方法和马氏链蒙特卡洛方法。以设备需求失效的稀少失效数据样本为例，推导了带Kass-Steffey修正的Beta-Binomial模型原理并编程求解，研究了马氏链蒙特卡罗方法及软件计算，对比了核电站后验失效概率的计算结果。计算表明：两种方法得到的部分失效概率后验估计的均值相差0～25%；95分位值相差5%～15%，两种方法都可用于稀少数据的样本估计。

参数经验贝叶斯方法； Kass-Steffey修正；MCMC；WinBUGS

设备可靠性参数是核电厂设备可靠性管理的重要内容，其在设备管理、老化管理、备品备件、维修规则等领域广泛应用。目前国内对于如何从多个核电站统计的失效数据直接得到各特定核电站设备可靠性参数的计算方法研究尚少[1]。随着我国《中国核电厂设备可靠性数据报告》(2015年)的发布[2]，已经具备采集并统计我国各核电站设备失效的数据样本的能力，亟须进一步完善有关利用统计的失效数据计算可靠性参数的贝叶斯方法，以得到我国的特定电站设备的可靠性参数。

本文研究了可靠性参数估计通常采用的分层模型原理[3，4]，以及实现分层模型的两种方法：带Kass-Steffey修正的参数经验贝叶斯方法和马氏链蒙特卡罗方法。这两种方法都可以给出核行业级的可靠性参数，也可以计算电站级的可靠性参数。鉴于核电站的失效数据多是稀少失效数据样本，本文以设备需求失效的稀少失效数据样本参数估计为例，详细推导了带Kass-Steffey 修正的Beta-Binomial模型的数学原理并编程计算，同时研究了马氏链蒙特卡罗方法及软件计算，最后对比了两种方法计算的特定电厂失效概率的后验估计值。

1 分层模型(Hierarchical Model)

m个核电站的某类失效数据若在部件类型、边界定义、运行状态、运行环境等方面相似，则失效数据xi可组成样本x=(x1，…，xm)，各个核电站失效数据的待估参数通常采用2层的分层模型进行估计[4]。

第1层：由于失效数据的相似性，将m个核电站的失效数据看成一个整体。尽管各个电厂的待估参数不同，但各个待估参数服从总体分布g。通过调用g分布m次，独立产生各个核电站参数。将g分布作为先验分布，g分布中的参数λ称为超参数。

第2层：依据各个核电站参数及失效数据分布，产生各个核电站的失效数据xi。

分层模型通过反复迭代，最后给出超参数λ和待估参数的后验估计值。

2 带Kass-Steffey修正的Beta-Binomial模型

Beta-Binomial模型用于描述设备的需求失效，例如阀门的开启/关闭失效、泵的启动失效、开关的打开/关闭失效、断路器的断开失效等。设第i个核电站某设备的失效次数为xi，需求次数为ni，失效概率pi，i=1, …,m，由于设备在每次需求时只有成功或失效2个结果，则xi服从参数pi的二项分布binomial(ni，pi)。

带Kass-Steffey修正的参数经验贝叶斯方法中，第一层：首先依据m个核电站的失效数据xi，采用最大似然法估计g分布的超参数。第二层：将g分布作为先验分布，使用贝叶斯方法，结合单个电厂数据，得到特定电厂参数pi的后验估计值[5, 6]。

2.1 Beta-Binomial模型的超参数求解

鉴于失效概率pi服从二项分布的共轭分布，因此总体分布g(p|α,β)服从beta(α,β)分布，即g的概率密度函数为：

(1)

因此，x的无条件分布服从Beta-Binomial分布，即：

(2)

m个核电站的似然函数为xi的无条件分布函数之积：

(3)

因此，似然函数的对数为：

(4)

其中，(4)式利用了伽马函数性质Γ(a)=(a-1)Γ(a-1)。为了便于似然函数的求导计算，需进行参数变换如下：

(5)

则(4)式变换为：

(6)

似然函数对数的偏导数方程组为：

(7)

超参数的估计值通过求解似然函数对数的偏导数方程组得到，即：

(8)

由于(8)式没有解析解，通常的算法[7]：先设置δ初始值，利用False position求根算法，求解方程组第一式中的μ值。然后将μ值带入方程组第二式，利用二分法再调整δ值。调整后的δ值再次作为初始值，带入第一式求解μ。如此反复迭代，直到方程组都收敛，最后给出μ和δ。

2.2 Beta-Binomial模型的特定电站参数pi估计

一阶修正下特定电厂参数pi的后验期望值不变[4, 5]，因此pi的后验期望值为：

(9)

带Kass-Steffey 修正的pi后验方差估计为：

(10)

其中，

(11)

同时，由于二阶协方差矩阵等于Fisher矩阵的逆，则:

var(μ)=1/J11

var(α)=1/J22

cov(μ,δ)=-J12/D

(12)

Fisher矩阵元为：

(13)

可以看到，带Kass-Steffey修正的参数经验贝叶斯方法是通过最大似然函数求导数解得超参数，然后利用超参数和特定电厂失效样本xi，计算特定电厂参数pi。值得注意的是，带Kass-Steffey修正的参数经验贝叶斯方法中的总体分布g选取的是二项分布的共轭分布。

3 马氏链蒙特卡罗方法

3.1 分层贝叶斯方法及MCMC

分层贝叶斯方法是完全的贝叶斯方法，不同于先计算超参数，对于需求失效数据，分层贝叶斯方法是将2个超参数和各个电厂待估参数组成维数为2+m维的先验参数矢量θ。θ的先验分布先假定：包括超参数α、β的超先验分布，m个核电站的pi条件分布及总体分布g，g分布可以假设为共轭分布，也可以是非共轭分布。

马氏链蒙特卡罗方法(MarkovChainMonteCarlo，MCMC)提供了一种计算分层贝叶斯模型的方法，其利用计算机模拟从后验分布抽样然后用样本估计值推断参数。Gibbs抽样是一种特殊的马氏链蒙特卡罗算法，具体是从某个初始点出发，通过全条件分布的循环抽样产生马氏链，算法如下[4]：

3.2 BUGS的Doodle模型的建立

BUGS(Bayesian inference Using Gibbs Sampling)是采用MCMC算法(利用Gibbs抽样)解决贝叶斯估计的软件。在WinBUGS软件中根据(1)式，在Doodle中建立相应的节点模型，如图1所示[8]。

图1 需求失效的Doodle模型Fig.1 The Doodle model of failures on demand

4 实例计算及结果比较

本文以68个核电站辅助给水系统电动部分启动失效的参数计算为例[4](累计失效个数为6个)，样本数据见表1，分别用Kass-Steffey 修正的Beta-Binomial模型和WinBUGS的方法进行计算，比较计算结果及比较见表2。

表1 AFW系统电动部分启动失效数据

表2 部分核电站后验需求失效概率的两种方法计算结果比较

在WinBUGS的Gibbs抽样模型中，假设超参数α服从均值和标准差都为1的指数超先验分布，β服从gamma(1.0, 0.035)超先验分布，p(i)为第i个核电站辅助给水系统电动部分启动失效概率，n(i)描述需求次数，x(i)描述失效次数，i=1,, 68。在文本中输入失效数据及需求次数数据，编译后加载超参数的初始值。在Sample Monitor Tool的节点中设置监测变量alpha、beta及p[i]，在Update中设置迭代次数。迭代完成后，在node中输入alpha、beta，查看stats获得参数估计结果。计算先验分布超参数计算的部分核电站后验需求失效概率见表2第3列。需要注意的是，WinBUGS软件使用时需要判断马氏链的收敛性、设置合理的迭代次数以及合适的超参数的超先验分布。

选取具有代表性的6个核电站(有失效数据次数的5个核电站和1个有0失效次数的核电站)，比较两种方法得到的后验估计值相对变化量见表2第4列。可以看到，两种方法得到的失效概率后验估计的均值和95分位值相差较小，均值相差0～25%，95分位值相差5%～15%；5分位值相差大，达32%～348%。

可以看到，对于后验参数的分位值估计，分位值越低，对模型假设的依赖越强，这一结论对于稀少样本的参数估计尤为明显，不同模型的计算结果差异越大。然而，在风险评价中，更多关注均值和95分位值(均值常作为点估计值，95分位值常作为风险估计的上限)，不同模型计算的低风险值(5分位值)的差异一定程度上可以不考虑[4]，这也是通常的可靠性参数采用均值和误差因子，而误差因子(Error Factor)通常采用95分位值除以均值计算的原因。因此，带Kass-Steffey修正的参数经验贝叶斯方法，也通常可以作为完全贝叶斯方法——马氏链蒙特卡罗方法的一阶近似，广泛应用于美国NRC的核电站的可靠性参数估计中[3]。此外，采用带Kass-Steffey修正的参数经验贝叶斯方法可以避免马氏链蒙特卡罗方法在收敛性判断、迭代次数、超先验分布设置合理性等方面的问题。

5 结论

本文以设备需求失效的稀少失效数据样本为例，分别用带Kass-Steffey修正的参数经验贝叶斯方法和马氏链蒙特卡罗方法计算了特定核电厂失效概率的后验估计值。计算表明：两种方法得到的部分核电厂失效概率后验估计的均值相差0～25%，95分位值相差5%～15%。对于大量失效数据的可靠性参数计算，两种方法的计算结果更为接近。一定程度上，带Kass-Steffey修正的参数经验贝叶斯方法可以作为分层贝叶斯马氏链蒙特卡洛方法的一阶近似。需要注意的是，带Kass-Steffey 修正的参数经验贝叶斯方法由于没有解析解，需要进行数值方法求解，有时可能因为数据的特殊性，存在无法求解的可能，且总体分布选取二项分布的共轭分布。马氏链蒙特卡罗方法在可靠性参数计算时则需要讨论收敛性判断、超先验分布设置的合理性、迭代次数等问题。这些问题都需要进一步深入研究。

[1] 茆定远, 薛大知.核电站PSA分析中可靠性数据处理的贝叶斯方法[J].核动力工程, 2000, 21(5):451-455.[2] 国家核安全局.中国核电厂设备可靠性数据报告，2015.

[3] Eide S A, Wierman T E, Gentillon C D, et al.Industry-Average Performance for Components and Initiating Events at U.S.Commercial Nuclear Power Plants [R].NUREG/CR-6928.2007.

[4] Atwood C L, LaChance J L, Martz H F, et al.Handbook of Parameter Estimation for Probabilistic Risk Assessment[R].NUREG/CR-6823.2003.

[5] Robert E.Kass, Duane Steffey, Approximate Bayesian Inference in Conditionally Independent Hierarchical Models (Parametric Empircal Bayes Models) [J].Journal of the American Statistical Association, 1989, 84(407)717-726.

[6] 陈妍, 郑鹏, 李朝君, 等.PSA 分析中可靠性参数的Kass-Steffey修正原理及应用研究[J].核动力工程, 2015, 36(3):70-74.

[7] Atwood, C.L. .1994, Hits per Trial: Basic Analysis of Binomial Data, EGG-RAAM-1 1041, Idaho National Engineering and Environmental Laboratory, Idaho Falls, ID.

[8] 陈妍, 余少青.WinBUGS软件在核电站可靠性参数估计中的应用[J], 科学技术与工程.2016, 16(1):151-154.

Calculation of Reliability Parameters of Nuclear Power Plant with Hierarchical Model

CHEN Yan，HE Liang，YU Shao-qing，CHEN Hai-ying，ZHANG Chun-ming

(Nuclear and Radiation Safety Center, Ministry of Environmental Protection, Beijing 100082)

There is little research on the calculation methods of evaluation reliability parameters directly from our own statistics failure data.The parametric empirical Bayes model with Kass-Steffey adjustment and the MCMC method are two Bayesian methods to estimate the reliability parameters.Take the failures on demand with sparse data as an example, this paper derives the principle of Beta-Binomial with Kass-Steffey adjustment, adopts two methods to calculate the hyperparameter and plant-specific posterior parameter.It shows that posterior mean of NPPs change 0～25%, posterior 95thchange 5%～15%.Two methods are all suitable for estimating the reliability parameters from sparse failure data.

Parametric empirical Bayes models; Kass-Steffey adjustment; MCMC; WinBUGS

2017-04-11

国家科技重大专项(2013ZX06002001-008)项目资助

陈 妍(1982—)，女，江苏徐州人，高级工程师，博士，理论物理，现从事核电站概率安全评价方面研究

何 亮：heliang@chinansc.cn

TL364

A

0258-0918(2017)03-0458-06