基于貝葉斯更新方法的235U熱中子獨立裂變產(chǎn)額協(xié)方差估計

王毅箴，崔夢蕾，郭 炯，鄔穎杰，劉保坤，李 富

(清華大學(xué) 核能與新能源技術(shù)研究院 先進核能技術(shù)協(xié)同創(chuàng)新中心先進反應(yīng)堆工程與安全教育部重點實驗室，北京100084)

裂變產(chǎn)額描述了裂變產(chǎn)物在裂變系統(tǒng)發(fā)生裂變反應(yīng)過程中產(chǎn)生的份額。根據(jù)裂變反應(yīng)的進程不同，裂變產(chǎn)額可分為獨立裂變產(chǎn)額和累積裂變產(chǎn)額。獨立裂變產(chǎn)額描述了瞬發(fā)裂變中子釋放后，在裂變產(chǎn)物發(fā)生衰變前各個產(chǎn)物產(chǎn)生的份額；累積裂變產(chǎn)額則描述了裂變產(chǎn)物發(fā)生長時間衰變后裂變產(chǎn)物的份額[1]。由于裂變過程的時間尺度小、復(fù)雜性高，所以裂變產(chǎn)額的測量與模擬具有較大的不確定度。目前，ENDF/B-VII.1數(shù)據(jù)庫中的235U熱中子裂變產(chǎn)額數(shù)據(jù)來自于1994年England等的評價結(jié)果[2]，該結(jié)果僅提供裂變產(chǎn)額各自的方差信息，不包含裂變產(chǎn)額間的協(xié)方差信息，導(dǎo)致采用基于抽樣統(tǒng)計不確定度分析方法研究裂變產(chǎn)額的不確定度傳遞時，無法給出合理自洽的裂變產(chǎn)額概率分布[3]。因此，裂變產(chǎn)額間的協(xié)方差估計與評價工作，十分必要[4-5]。

貝葉斯更新方法廣泛應(yīng)用于數(shù)據(jù)同化、數(shù)據(jù)調(diào)整及模型擬合等領(lǐng)域。在對裂變產(chǎn)額協(xié)方差的研究中， Kawano等首先基于鏈裂變產(chǎn)額估計了239Pu裂變產(chǎn)額的協(xié)方差[6]。其次， Pigni等將Kawano等的方法應(yīng)用于ENDF/B-VII.1數(shù)據(jù)庫中，研究了基于累積裂變產(chǎn)額的獨立裂變產(chǎn)額調(diào)整方法，但未重點研究更新后的獨立裂變產(chǎn)額的協(xié)方差矩陣特征[7]。隨后， Fiorito等研究評價了分別用鏈裂變產(chǎn)額和累積裂變產(chǎn)額估計獨立裂變產(chǎn)額協(xié)方差的差異及其對傳遞的裂變產(chǎn)額不確定度的影響[8-9]。本文基于貝葉斯更新方法，利用ENDF/B-VII.1數(shù)據(jù)庫中的累積裂變產(chǎn)額數(shù)據(jù)及獨立裂變產(chǎn)額自身的物理約束條件，對ENDF/B-VII.1數(shù)據(jù)庫中獨立裂變產(chǎn)額的協(xié)方差進行估計，并對估計結(jié)果進行分析和討論。

1 獨立裂變產(chǎn)額的協(xié)方差估計方法

1.1 貝葉斯更新方法

貝葉斯更新方法(Bayesian updating method)是一種核數(shù)據(jù)調(diào)整方法，它根據(jù)核數(shù)據(jù)的新認知(新評估結(jié)果或物理約束條件)對核數(shù)據(jù)進行調(diào)整，從而保持核數(shù)據(jù)與新認知的一致性[10-11]。該調(diào)整過程是在概率空間定義下進行的，因此，核數(shù)據(jù)x先被賦予某種概率密度分布，稱為先驗概率密度函數(shù)p(x)。核數(shù)據(jù)的新認知信息y被同樣賦予某種概率密度分布，從而計算似然值p(y|x)。根據(jù)貝葉斯定理，核數(shù)據(jù)的后驗概率分布p(x|y)可表示為

(1)

后驗概率分布的最大值即為核數(shù)據(jù)的最佳調(diào)整值。由式(1)可知，貝葉斯更新方法對核數(shù)據(jù)進行調(diào)整具有3個特點：

1)允許按序分級地引入不同的核數(shù)據(jù)新認知y1,y2,…,yN，這樣可以避免每當有新認知產(chǎn)生時，對以往核數(shù)據(jù)認知的重新引入與計算，也可以方便地分級觀察和研究不同新認知對核數(shù)據(jù)調(diào)整值的影響。式(1)可改寫為

(2)

2)貝葉斯更新給出的核數(shù)據(jù)最佳調(diào)整值依賴于先驗概率密度分布和似然函數(shù)，它們可以通過以往人們對核數(shù)據(jù)概率分布的研究結(jié)果或基于最大信息熵原理給出[10]，從而保證核數(shù)據(jù)調(diào)整的合理性。

3)由于貝葉斯更新方法是基于核數(shù)據(jù)后驗概率分布進行調(diào)整的，因此，基于獲得的后驗概率分布，可以對核數(shù)據(jù)的高階矩信息進行模擬或直接計算，從而可以獲得核數(shù)據(jù)在新認知約束下的協(xié)方差信息。

本文基于貝葉斯更新方法對獨立裂變產(chǎn)額進行調(diào)整，并對ENDF/B-VII.1核數(shù)據(jù)庫中缺失的協(xié)方差信息進行估計。由于ENDF/B-VII.1核數(shù)據(jù)庫中僅給出了獨立裂變產(chǎn)額的均值和標準差，由最大信息熵原理可知，采用正態(tài)分布描述獨立裂變產(chǎn)額能夠得到最大信息熵，從而可以充分利用已知的不確定度信息，最小地引入人為的主觀影響。同理，對累積裂變產(chǎn)額也同樣賦予正態(tài)分布計算似然值。

(3)

(4)

其中，Yi為獨立裂變產(chǎn)額向量，Yi∈RP×1；P為總裂變產(chǎn)物數(shù)目；Yc,0,Yi,0分別為ENDF/B-VII.1提供的累積裂變產(chǎn)額與獨立裂變產(chǎn)額的均值；Vc,Vi,0分別為累積裂變產(chǎn)額和獨立裂變產(chǎn)額的協(xié)方差矩陣，它們均為對角陣，對角元素為ENDF/B-VII.1數(shù)據(jù)庫中提供的各個裂變產(chǎn)額的方差,Vc∈RP×P，Vi,0∈RP×P；M為累積裂變產(chǎn)額與獨立裂變產(chǎn)額的線性映射矩陣。每個裂變產(chǎn)物核素可由其質(zhì)量數(shù)A、質(zhì)子數(shù)Z及其同核異能態(tài)I確定，記為(A,Z,I)。累積裂變產(chǎn)額與獨立裂變產(chǎn)額之間的關(guān)系可表示為[1]

Yc(A,Z,I)=Yi(A,Z,I)+

(5)

其中，Yc(A,Z,I),Yi(A,Z,I)分別表示裂變產(chǎn)物核素(A,Z,I)的累積裂變產(chǎn)額和獨立裂變產(chǎn)額；b(A′,Z′,I′→A,Z,I)為裂變產(chǎn)物核素(A′,Z′,I′)發(fā)生衰變到裂變產(chǎn)物(A,Z,I)的分支比。聯(lián)立所有裂變產(chǎn)物核素的累積裂變產(chǎn)額計算式，根據(jù)式(5)可得

Yc=Yi+BYc

(6)

Yc=(E-B)-1Yi

(7)

Yc=MYi

(8)

其中，Yc為累積裂變產(chǎn)額向量，Yc∈RP×1；B為裂變產(chǎn)物核素之間的分支比矩陣，B∈RP×P；E為單位矩陣，E∈RP×P。

由式(1)可以得到累積裂變產(chǎn)額更新后的獨立裂變產(chǎn)額的后驗概率分布為

(9)

經(jīng)推導(dǎo)可得，獨立裂變產(chǎn)額的后驗概率分布也為正態(tài)分布，只是分布參數(shù)發(fā)生了變化。因此，獨立裂變產(chǎn)額的后驗概率分布可表示為

(10)

其中，

(11)

Yi,1=Yi,0+Vi,1MTV-1(Yc,0-MYi,0)

(12)

由于正態(tài)分布均值點的鄰域內(nèi)具有最大概率，所以，獨立裂變產(chǎn)額的最佳調(diào)整值即為后驗概率分布的均值Yi,1。后驗概率分布的協(xié)方差矩陣Vi,1，可用于研究獨立裂變產(chǎn)額不確定度的傳遞。

1.2 獨立裂變產(chǎn)額的物理約束條件

獨立裂變產(chǎn)額應(yīng)滿足3個物理約束條件：

1)在忽略三元裂變的情況下，總獨立裂變產(chǎn)額應(yīng)為2.0；

2)在忽略三元裂變的情況下，各個裂變產(chǎn)物的質(zhì)量數(shù)以其各自的獨立裂變產(chǎn)額為權(quán)重的平均值，應(yīng)為裂變系統(tǒng)裂變前復(fù)合核的總質(zhì)量減去釋放的平均瞬發(fā)中子數(shù)；

3)各個裂變產(chǎn)物的電荷數(shù)以其各自的獨立裂變產(chǎn)額為權(quán)重的平均值，應(yīng)與裂變系統(tǒng)的總電荷數(shù)相等。

由于本文依據(jù)的ENDF/B-VII.1數(shù)據(jù)庫中未提供H，He, Li, Be, B, C, N, Ne, Ti等輕核素的相關(guān)數(shù)據(jù)[12]，所以，本文分析中均忽略三元裂變。

總裂變產(chǎn)額維持二裂變的約束條件，采用正態(tài)分布計算似然值，則

p(Yi|Ytot,Yc,0)∝p(Ytot|Yi,Yc,0)p(Yi|Yc,0)∝

(13)

經(jīng)推導(dǎo)，引入二裂變守恒后的獨立裂變產(chǎn)額后驗概率分布p(Yi|Ytot,Yc,0)，可化為正態(tài)分布形式，即

(14)

其中，

(15)

(16)

其中，Ytot為總獨立裂變產(chǎn)額；σ2為總裂變產(chǎn)額的求和精度，本文取σ2=10-5；U為單位向量，U∈RP×1。

同理，引入裂變系統(tǒng)總質(zhì)量數(shù)守恒約束條件，則

p(Yi│Atot,Ytot,Yc,0)∝p(Atot│Yi,Ytot,Yc,0)·

(17)

p(Yi│Atot,Ytot,Yc,0)∝

(18)

(19)

(20)

其中，Atot為裂變系統(tǒng)總質(zhì)量數(shù)，該值可通過各個裂變核素的質(zhì)量數(shù)以其ENDF/B-VII.1數(shù)據(jù)庫中提供的獨立裂變產(chǎn)額加權(quán)平均給出；A為各個裂變核素的質(zhì)量數(shù)向量，A∈RP×1。

引入裂變系統(tǒng)的總電荷數(shù)守恒約束條件，則

p(Yi│Ztot,Atot,Ytot,Yc,0)∝p(Ztot│Yi,Atot,Ytot,Yc,0)·

(21)

p(Yi│Ztot,Atot,Ytot,Yc,0)∝

(22)

(23)

(24)

其中，Ztot為裂變系統(tǒng)總電荷數(shù)，該值可通過各個裂變核素的電荷數(shù)以其ENDF/B-VII.1數(shù)據(jù)庫中提供的獨立裂變裂變產(chǎn)額加權(quán)平均給出；Z為各個裂變核素的電荷數(shù)向量，Z∈RP×1。

最終得到的獨立裂變產(chǎn)額的后驗概率分布p(Yi│Ztot,Atot,Ytot,Yc,0)，能夠給出獨立裂變產(chǎn)額在累積裂變產(chǎn)額更新及自身的物理約束下的最佳調(diào)整值Yi,4及合理估計的協(xié)方差矩陣Vi,4。利用Yi,4和Vi,4進行獨立裂變產(chǎn)額抽樣，能夠獲得符合獨立裂變物理約束條件且累積裂變產(chǎn)額一致的獨立裂變產(chǎn)額樣本數(shù)據(jù)，為后續(xù)開展不確定度傳遞研究提供基礎(chǔ)。

2 估計結(jié)果

2.1 獨立裂變產(chǎn)額調(diào)整結(jié)果與討論

獨立裂變產(chǎn)額的調(diào)整結(jié)果，如圖1所示。其中，藍色點代表ENDF/B-VII.1數(shù)據(jù)庫中提供的235U熱中子裂變的獨立裂變產(chǎn)額，即先驗獨立裂變產(chǎn)額；橙色點代表經(jīng)過貝葉斯更新后的獨立裂變產(chǎn)額最佳調(diào)整值Yi,4，即后驗獨立裂變產(chǎn)額的均值。為方便比較分析，縱軸給出的是獨立裂變質(zhì)量產(chǎn)額Yi,m，即按質(zhì)量數(shù)對裂變核素分組，并對每個質(zhì)量組內(nèi)核素的獨立裂變產(chǎn)額求和的結(jié)果。由圖1可見，獨立裂變質(zhì)量產(chǎn)額調(diào)整較大的核素在曲線的兩峰處，對應(yīng)質(zhì)量數(shù)分別為80～110和130～150。這是由于兩峰處裂變產(chǎn)物累積裂變產(chǎn)額的測量不確定度較小，累積裂變產(chǎn)額較其他位置核素的累積裂變產(chǎn)額更為準確，通過貝葉斯更新方法對兩峰處裂變產(chǎn)物的獨立裂變產(chǎn)額調(diào)整時，調(diào)整效果較為明顯。

圖1 獨立裂變產(chǎn)額的調(diào)整結(jié)果Fig.1 Independent fission yields adjustment results

為了分析貝葉斯更新調(diào)整后的獨立裂變產(chǎn)額Yi,4與ENDF/B-VII.1數(shù)據(jù)庫提供的累積裂變產(chǎn)額評價值Yc，0的一致性，根據(jù)式(8)，分別利用先驗獨立裂變產(chǎn)額Yi，0與后驗獨立裂變產(chǎn)額Yi,4計算對應(yīng)的累積裂變產(chǎn)額，并記為先驗累積裂變產(chǎn)額計算值與后驗累積裂變產(chǎn)額計算值。為方便分析，將累積裂變產(chǎn)額按照質(zhì)量數(shù)進行分組求和，得到先驗累積裂變質(zhì)量產(chǎn)額與后驗累積裂變質(zhì)量產(chǎn)額計算值。然后分別計算它們與ENDF/B-VII.1數(shù)據(jù)庫中提供的累積裂變產(chǎn)額評價值之比r(Yc,m)，得到r(Yc,m)隨質(zhì)量數(shù)A的分布曲線，如圖2所示。

由圖2可見，用后驗獨立裂變產(chǎn)額計算的累積裂變產(chǎn)額的估計不確定度在20%以內(nèi)，與用先驗獨立裂變產(chǎn)額計算的累積裂變產(chǎn)額的估計不確定度相比減少約30%，提高了ENDF/B-VII.1獨立裂變產(chǎn)額與累積裂變產(chǎn)額的一致性。對質(zhì)量數(shù)為110～130的核素，更新后的估計不確定度有所放大，這是由于這些核素的累積裂變產(chǎn)額的先驗計算結(jié)果與ENDF/B-VII.1數(shù)據(jù)庫提供的評價值之間的偏差較大造成的。為表征先驗累積裂變產(chǎn)額計算值與ENDF/B-VII.1數(shù)據(jù)庫提供的累積裂變產(chǎn)額評價值的一致性，定義一致性因子H為

(25)

其中，NA為質(zhì)量鏈A中裂變產(chǎn)物核素的數(shù)目；Yc,n為質(zhì)量鏈A中裂變產(chǎn)物n先驗累積裂變產(chǎn)額計算值；Yc,0,n為質(zhì)量鏈A中ENDF/B-VII.1提供的裂變產(chǎn)物n累積裂變產(chǎn)額評價值；σc,0,n為ENDF/B-VII.1提供的裂變產(chǎn)物n累積裂變產(chǎn)額的標準偏差[13]。H越大，表示計算值與評價值之間的一致性越弱。圖3對比了先驗累積裂變質(zhì)量產(chǎn)額計算值與ENDF/B-VII.1數(shù)據(jù)庫中評價值的一致性。由圖3可見，質(zhì)量數(shù)為110～130和80～90核素的一致性因子H較其周圍核素的H值大，具有較弱的一致性，這種較弱的一致性影響了貝葉斯更新的獨立裂變產(chǎn)額調(diào)整結(jié)果。由于這些核素的獨立裂變產(chǎn)額較小(位于產(chǎn)額分布曲線上兩峰之間的低谷處)，因此，這些核素的獨立裂變產(chǎn)額可能需要在未來的評價數(shù)據(jù)庫中重新評估。

圖3 先驗累積裂變質(zhì)量產(chǎn)額計算值與ENDF/B-VII.1數(shù)據(jù)庫中評價值的一致性對比Fig.3 Prior calculated cumulative fission yields consistencycomparison with ENDF/B-VII.1 evaluations

圖4給出了后驗與先驗獨立裂變產(chǎn)額標準偏差的比值r(σ)隨質(zhì)量數(shù)A及電荷數(shù)Z的分布。由圖4可見，多數(shù)裂變產(chǎn)物獨立裂變產(chǎn)額的不確定度經(jīng)過更新后有所減小，約為更新前的70%，而且位于兩峰處的核素，獨立裂變產(chǎn)額不確定度減小的幅度較大，約為更新前的20%。這是由于隨著累積裂變產(chǎn)額評價值及其引入的不確定度的減小，且經(jīng)過獨立裂變產(chǎn)額物理約束條件的限制和似然函數(shù)的加權(quán)調(diào)整，獨立裂變產(chǎn)額的后驗概率分布會在最佳調(diào)整值附近聚集更多，從而使后驗概率分布函數(shù)更加尖銳，標準偏差更小。

圖4 后驗與先驗獨立裂變產(chǎn)額標準偏差的比值r(σ)隨質(zhì)量數(shù)A及電荷數(shù)Z的分布Fig.4 Independent fission yields posterior to priorstandard deviation ratios r(σ) distribution againstmass number A and charge number Z

2.2 獨立裂變產(chǎn)額相關(guān)性估計結(jié)果與討論

隨著累積裂變產(chǎn)額測量值及獨立裂變產(chǎn)額物理約束條件的引入，獨立裂變產(chǎn)額后驗協(xié)方差矩陣會不斷發(fā)生變化。圖5給出了不同更新步驟下獨立裂變產(chǎn)額相關(guān)系數(shù)矩陣的變化過程。本文對裂變產(chǎn)物按照質(zhì)量數(shù)大小進行分組并降序排列，在每個質(zhì)量數(shù)分組內(nèi)按照裂變產(chǎn)物的電荷數(shù)大小升序排列。

由圖5可見，相關(guān)系數(shù)呈團簇狀出現(xiàn)在主對角線附近，這是由于累積裂變產(chǎn)額和獨立裂變產(chǎn)額之間是通過裂變產(chǎn)物的放射性衰變相聯(lián)系的，相關(guān)性主要出現(xiàn)在有放射性衰變聯(lián)系的核素之間。隨著二裂變約束的引入，可以看到圖中Vi,2在相距較遠的核素位置處產(chǎn)生了相關(guān)性。而隨著裂變系統(tǒng)的質(zhì)量數(shù)和電荷數(shù)守恒的引入，更多的相關(guān)性被引入到獨立裂變產(chǎn)額的協(xié)方差矩陣中。

(a)Vi,1:p(Yi|Yc,0)

(b)Vi,2:p(Yi|Ytot,Yc,0)

(c)Vi,3:p(Yi|Atot,Ytot,Yc,0)

(d)Vi,4:p(Yi|Ztot,Atot,Ytot,Yc,0)

圖6為圖5中的Vi,4的局部放大圖。其中，紅色點代表正相關(guān)性系數(shù)，黑色點代表負相關(guān)性系數(shù)。由圖6可見，在主對角線附近呈團簇狀出現(xiàn)的多為負相關(guān)性系數(shù)，即由裂變產(chǎn)物衰變引起的核素間的相關(guān)性多為負相關(guān)；而在主對角線以外由獨立裂變產(chǎn)額的物理約束引起的相關(guān)性則多為正相關(guān)。圖7為圖6左側(cè)矩陣的下三角元素的相關(guān)系數(shù)散點圖。由圖7可見，貝葉斯更新方法得到的獨立裂變產(chǎn)額之間的相關(guān)性系數(shù)多為負相關(guān)，且多數(shù)的相關(guān)性較弱，相關(guān)系數(shù)小于0.25。少數(shù)核素之間存在較強的負相關(guān)性，極少數(shù)核素之間存在較強的正相關(guān)性。

圖6 相關(guān)系數(shù)矩陣Vi,4的局部放大圖Fig.6 Correlation matrix of independentfission yields Vi,4 with local zoom-in view

圖7 相關(guān)系數(shù)矩陣Vi,4 下三角元素散點圖Fig.7 Scatter plot of lower triangle elements in correlation matrix Vi,4

近年來，國外雖然已經(jīng)開展了很多基于實驗鏈產(chǎn)額和累積裂變產(chǎn)額數(shù)據(jù)更新獨立裂變產(chǎn)額和協(xié)方差數(shù)據(jù)的研究工作，但是，目前ENDF/B-VII.1數(shù)據(jù)庫還未提供獨立裂變產(chǎn)額的協(xié)方差數(shù)據(jù)。文獻[14]中利用本文研究的獨立裂變產(chǎn)額的協(xié)方差數(shù)據(jù)，研究了球床式高溫氣冷堆裂變產(chǎn)額的不確定度傳遞及獨立裂變產(chǎn)額不確定度對堆芯平衡態(tài)下有效增殖因子和核素濃度不確定度的貢獻，結(jié)果表明，有效增殖因子的相對不確定度約為2.44×10-4(燃耗為90 GW·d·t-1)，該結(jié)果和REBUS基準題中單一燃料棒柵元的相應(yīng)評價結(jié)果6.0×10-4(燃耗為54 GW·d·t-1)較為接近[9]，說明本文的評估工作是合理的。

3 結(jié)論

本文基于貝葉斯更新方法，分級逐步將ENDF/B-VII.1數(shù)據(jù)庫提供的累積裂變產(chǎn)額測量數(shù)據(jù)及獨立裂變產(chǎn)額服從的二裂變守恒、裂變系統(tǒng)質(zhì)量數(shù)守恒和電荷數(shù)守恒等條件引入?yún)f(xié)方差估計過程中，對ENDF/B-VII.1數(shù)據(jù)庫中235U熱中子裂變的獨立裂變產(chǎn)額協(xié)方差進行了估計。結(jié)果表明，與貝葉斯更新前相比，更新后多數(shù)核素獨立裂變產(chǎn)額的不確定度均出現(xiàn)了不同程度的減小，不確定度約減少了30%。多數(shù)核素的獨立裂變產(chǎn)額之間為負相關(guān)，且相關(guān)性較弱，極少數(shù)核素的獨立裂變產(chǎn)額之間存在較強的正相關(guān)。

本文研究是進一步開展獨立裂變產(chǎn)額不確定度傳遞研究工作的基礎(chǔ)，未來將基于本文估計的獨立裂變產(chǎn)額協(xié)方差矩陣，研究裂變產(chǎn)額不確定度對反應(yīng)堆燃耗及衰變熱不確定度的貢獻。