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

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

(環境保護部核與輻射安全中心，北京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