陳 妍，何 亮，余少青，陳海英，張春明

(環(huán)境保護部核與輻射安全中心，北京100082)

?

應用分層模型進行核電站設備可靠性參數(shù)估計

陳 妍，何 亮，余少青，陳海英，張春明

(環(huán)境保護部核與輻射安全中心，北京100082)

目前國內對于從多個核電站統(tǒng)計的設備失效數(shù)據(jù)進行各特定核電站設備可靠性參數(shù)估計的方法研究尚少。本文研究了用于可靠性參數(shù)估計的分層模型以及實現(xiàn)分層模型的兩種方法：帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法和馬氏鏈蒙特卡洛方法。以設備需求失效的稀少失效數(shù)據(jù)樣本為例，推導了帶Kass-Steffey修正的Beta-Binomial模型原理并編程求解，研究了馬氏鏈蒙特卡羅方法及軟件計算，對比了核電站后驗失效概率的計算結果。計算表明：兩種方法得到的部分失效概率后驗估計的均值相差0～25%；95分位值相差5%～15%，兩種方法都可用于稀少數(shù)據(jù)的樣本估計。

參數(shù)經(jīng)驗貝葉斯方法； Kass-Steffey修正；MCMC；WinBUGS

設備可靠性參數(shù)是核電廠設備可靠性管理的重要內容，其在設備管理、老化管理、備品備件、維修規(guī)則等領域廣泛應用。目前國內對于如何從多個核電站統(tǒng)計的失效數(shù)據(jù)直接得到各特定核電站設備可靠性參數(shù)的計算方法研究尚少[1]。隨著我國《中國核電廠設備可靠性數(shù)據(jù)報告》(2015年)的發(fā)布[2]，已經(jīng)具備采集并統(tǒng)計我國各核電站設備失效的數(shù)據(jù)樣本的能力，亟須進一步完善有關利用統(tǒng)計的失效數(shù)據(jù)計算可靠性參數(shù)的貝葉斯方法，以得到我國的特定電站設備的可靠性參數(shù)。

本文研究了可靠性參數(shù)估計通常采用的分層模型原理[3，4]，以及實現(xiàn)分層模型的兩種方法：帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法和馬氏鏈蒙特卡羅方法。這兩種方法都可以給出核行業(yè)級的可靠性參數(shù)，也可以計算電站級的可靠性參數(shù)。鑒于核電站的失效數(shù)據(jù)多是稀少失效數(shù)據(jù)樣本，本文以設備需求失效的稀少失效數(shù)據(jù)樣本參數(shù)估計為例，詳細推導了帶Kass-Steffey 修正的Beta-Binomial模型的數(shù)學原理并編程計算，同時研究了馬氏鏈蒙特卡羅方法及軟件計算，最后對比了兩種方法計算的特定電廠失效概率的后驗估計值。

1 分層模型(Hierarchical Model)

m個核電站的某類失效數(shù)據(jù)若在部件類型、邊界定義、運行狀態(tài)、運行環(huán)境等方面相似，則失效數(shù)據(jù)xi可組成樣本x=(x1，…，xm)，各個核電站失效數(shù)據(jù)的待估參數(shù)通常采用2層的分層模型進行估計[4]。

第1層：由于失效數(shù)據(jù)的相似性，將m個核電站的失效數(shù)據(jù)看成一個整體。盡管各個電廠的待估參數(shù)不同，但各個待估參數(shù)服從總體分布g。通過調用g分布m次，獨立產(chǎn)生各個核電站參數(shù)。將g分布作為先驗分布，g分布中的參數(shù)λ稱為超參數(shù)。

第2層：依據(jù)各個核電站參數(shù)及失效數(shù)據(jù)分布，產(chǎn)生各個核電站的失效數(shù)據(jù)xi。

分層模型通過反復迭代，最后給出超參數(shù)λ和待估參數(shù)的后驗估計值。

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

Beta-Binomial模型用于描述設備的需求失效，例如閥門的開啟/關閉失效、泵的啟動失效、開關的打開/關閉失效、斷路器的斷開失效等。設第i個核電站某設備的失效次數(shù)為xi，需求次數(shù)為ni，失效概率pi，i=1, …,m，由于設備在每次需求時只有成功或失效2個結果，則xi服從參數(shù)pi的二項分布binomial(ni，pi)。

帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法中，第一層：首先依據(jù)m個核電站的失效數(shù)據(jù)xi，采用最大似然法估計g分布的超參數(shù)。第二層：將g分布作為先驗分布，使用貝葉斯方法，結合單個電廠數(shù)據(jù)，得到特定電廠參數(shù)pi的后驗估計值[5, 6]。

2.1 Beta-Binomial模型的超參數(shù)求解

鑒于失效概率pi服從二項分布的共軛分布，因此總體分布g(p|α,β)服從beta(α,β)分布，即g的概率密度函數(shù)為：

(1)

因此，x的無條件分布服從Beta-Binomial分布，即：

(2)

m個核電站的似然函數(shù)為xi的無條件分布函數(shù)之積：

(3)

因此，似然函數(shù)的對數(shù)為：

(4)

其中，(4)式利用了伽馬函數(shù)性質Γ(a)=(a-1)Γ(a-1)。為了便于似然函數(shù)的求導計算，需進行參數(shù)變換如下：

(5)

則(4)式變換為：

(6)

似然函數(shù)對數(shù)的偏導數(shù)方程組為：

(7)

超參數(shù)的估計值通過求解似然函數(shù)對數(shù)的偏導數(shù)方程組得到，即：

(8)

由于(8)式?jīng)]有解析解，通常的算法[7]：先設置δ初始值，利用False position求根算法，求解方程組第一式中的μ值。然后將μ值帶入方程組第二式，利用二分法再調整δ值。調整后的δ值再次作為初始值，帶入第一式求解μ。如此反復迭代，直到方程組都收斂，最后給出μ和δ。

2.2 Beta-Binomial模型的特定電站參數(shù)pi估計

一階修正下特定電廠參數(shù)pi的后驗期望值不變[4, 5]，因此pi的后驗期望值為：

(9)

帶Kass-Steffey 修正的pi后驗方差估計為：

(10)

其中，

(11)

同時，由于二階協(xié)方差矩陣等于Fisher矩陣的逆，則:

var(μ)=1/J11

var(α)=1/J22

cov(μ,δ)=-J12/D

(12)

Fisher矩陣元為：

(13)

可以看到，帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法是通過最大似然函數(shù)求導數(shù)解得超參數(shù)，然后利用超參數(shù)和特定電廠失效樣本xi，計算特定電廠參數(shù)pi。值得注意的是，帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法中的總體分布g選取的是二項分布的共軛分布。

3 馬氏鏈蒙特卡羅方法

3.1 分層貝葉斯方法及MCMC

分層貝葉斯方法是完全的貝葉斯方法，不同于先計算超參數(shù)，對于需求失效數(shù)據(jù)，分層貝葉斯方法是將2個超參數(shù)和各個電廠待估參數(shù)組成維數(shù)為2+m維的先驗參數(shù)矢量θ。θ的先驗分布先假定：包括超參數(shù)α、β的超先驗分布，m個核電站的pi條件分布及總體分布g，g分布可以假設為共軛分布，也可以是非共軛分布。

馬氏鏈蒙特卡羅方法(MarkovChainMonteCarlo，MCMC)提供了一種計算分層貝葉斯模型的方法，其利用計算機模擬從后驗分布抽樣然后用樣本估計值推斷參數(shù)。Gibbs抽樣是一種特殊的馬氏鏈蒙特卡羅算法，具體是從某個初始點出發(fā)，通過全條件分布的循環(huán)抽樣產(chǎn)生馬氏鏈，算法如下[4]：

3.2 BUGS的Doodle模型的建立

BUGS(Bayesian inference Using Gibbs Sampling)是采用MCMC算法(利用Gibbs抽樣)解決貝葉斯估計的軟件。在WinBUGS軟件中根據(jù)(1)式，在Doodle中建立相應的節(jié)點模型，如圖1所示[8]。

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

4 實例計算及結果比較

本文以68個核電站輔助給水系統(tǒng)電動部分啟動失效的參數(shù)計算為例[4](累計失效個數(shù)為6個)，樣本數(shù)據(jù)見表1，分別用Kass-Steffey 修正的Beta-Binomial模型和WinBUGS的方法進行計算，比較計算結果及比較見表2。

表1 AFW系統(tǒng)電動部分啟動失效數(shù)據(jù)

表2 部分核電站后驗需求失效概率的兩種方法計算結果比較

在WinBUGS的Gibbs抽樣模型中，假設超參數(shù)α服從均值和標準差都為1的指數(shù)超先驗分布，β服從gamma(1.0, 0.035)超先驗分布，p(i)為第i個核電站輔助給水系統(tǒng)電動部分啟動失效概率，n(i)描述需求次數(shù)，x(i)描述失效次數(shù)，i=1,, 68。在文本中輸入失效數(shù)據(jù)及需求次數(shù)數(shù)據(jù)，編譯后加載超參數(shù)的初始值。在Sample Monitor Tool的節(jié)點中設置監(jiān)測變量alpha、beta及p[i]，在Update中設置迭代次數(shù)。迭代完成后，在node中輸入alpha、beta，查看stats獲得參數(shù)估計結果。計算先驗分布超參數(shù)計算的部分核電站后驗需求失效概率見表2第3列。需要注意的是，WinBUGS軟件使用時需要判斷馬氏鏈的收斂性、設置合理的迭代次數(shù)以及合適的超參數(shù)的超先驗分布。

選取具有代表性的6個核電站(有失效數(shù)據(jù)次數(shù)的5個核電站和1個有0失效次數(shù)的核電站)，比較兩種方法得到的后驗估計值相對變化量見表2第4列?？梢钥吹剑瑑煞N方法得到的失效概率后驗估計的均值和95分位值相差較小，均值相差0～25%，95分位值相差5%～15%；5分位值相差大，達32%～348%。

可以看到，對于后驗參數(shù)的分位值估計，分位值越低，對模型假設的依賴越強，這一結論對于稀少樣本的參數(shù)估計尤為明顯，不同模型的計算結果差異越大。然而，在風險評價中，更多關注均值和95分位值(均值常作為點估計值，95分位值常作為風險估計的上限)，不同模型計算的低風險值(5分位值)的差異一定程度上可以不考慮[4]，這也是通常的可靠性參數(shù)采用均值和誤差因子，而誤差因子(Error Factor)通常采用95分位值除以均值計算的原因。因此，帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法，也通?？梢宰鳛橥耆惾~斯方法——馬氏鏈蒙特卡羅方法的一階近似，廣泛應用于美國NRC的核電站的可靠性參數(shù)估計中[3]。此外，采用帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法可以避免馬氏鏈蒙特卡羅方法在收斂性判斷、迭代次數(shù)、超先驗分布設置合理性等方面的問題。

5 結論

本文以設備需求失效的稀少失效數(shù)據(jù)樣本為例，分別用帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法和馬氏鏈蒙特卡羅方法計算了特定核電廠失效概率的后驗估計值。計算表明：兩種方法得到的部分核電廠失效概率后驗估計的均值相差0～25%，95分位值相差5%～15%。對于大量失效數(shù)據(jù)的可靠性參數(shù)計算，兩種方法的計算結果更為接近。一定程度上，帶Kass-Steffey修正的參數(shù)經(jīng)驗貝葉斯方法可以作為分層貝葉斯馬氏鏈蒙特卡洛方法的一階近似。需要注意的是，帶Kass-Steffey 修正的參數(shù)經(jīng)驗貝葉斯方法由于沒有解析解，需要進行數(shù)值方法求解，有時可能因為數(shù)據(jù)的特殊性，存在無法求解的可能，且總體分布選取二項分布的共軛分布。馬氏鏈蒙特卡羅方法在可靠性參數(shù)計算時則需要討論收斂性判斷、超先驗分布設置的合理性、迭代次數(shù)等問題。這些問題都需要進一步深入研究。

[1] 茆定遠, 薛大知.核電站PSA分析中可靠性數(shù)據(jù)處理的貝葉斯方法[J].核動力工程, 2000, 21(5):451-455.[2] 國家核安全局.中國核電廠設備可靠性數(shù)據(jù)報告，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 分析中可靠性參數(shù)的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軟件在核電站可靠性參數(shù)估計中的應用[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—)，女，江蘇徐州人，高級工程師，博士，理論物理，現(xiàn)從事核電站概率安全評價方面研究

何 亮：heliang@chinansc.cn

TL364

A

0258-0918(2017)03-0458-06