馬續(xù)波，李耀舟，吳 屈，彭星杰，黃自鋒，朱潤(rùn)澤，張 斌

(1.華北電力大學(xué) 核科學(xué)與工程學(xué)院，北京 102206；2.中國(guó)核動(dòng)力研究設(shè)計(jì)院 核反應(yīng)堆系統(tǒng)設(shè)計(jì)技術(shù)重點(diǎn)實(shí)驗(yàn)室，四川 成都 610213)

蒙特卡羅方法已廣泛應(yīng)用到了各領(lǐng)域。由于具有已知分布的隨機(jī)變量抽樣方法在蒙特卡羅方法中占有重要位置，并已得到廣泛研究[1-4]。采用統(tǒng)計(jì)抽樣方法進(jìn)行不確定度和靈敏度分析也成為了核工程領(lǐng)域研究熱點(diǎn)[5-7]，但采用統(tǒng)計(jì)抽樣法進(jìn)行不確定度分析時(shí)，需大量抽樣樣本，而對(duì)于有些問題，問題求解本身需消耗大量時(shí)間，如求解玻爾茲曼輸運(yùn)方程、熱傳導(dǎo)方程等，這使采用統(tǒng)計(jì)抽樣法進(jìn)行不確定度分析需消耗大量時(shí)間，限制了統(tǒng)計(jì)抽樣法的應(yīng)用和推廣。

采用統(tǒng)計(jì)抽樣法進(jìn)行不確定度分析，首先要產(chǎn)生多維隨機(jī)變量序列，而對(duì)于多維隨機(jī)變量的產(chǎn)生，需根據(jù)多維隨機(jī)變量之間的協(xié)方差矩陣，產(chǎn)生具有相關(guān)變量的隨機(jī)數(shù)序列。文獻(xiàn)[5]給出了相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法的理論推導(dǎo)及數(shù)值驗(yàn)證。雖從理論上講，產(chǎn)生的相關(guān)變量隨機(jī)數(shù)序列的協(xié)方差矩陣與真實(shí)的協(xié)方差矩陣應(yīng)一致，但研究發(fā)現(xiàn)，僅在大樣本情況下，產(chǎn)生的相關(guān)變量隨機(jī)數(shù)的協(xié)方差矩陣才與真實(shí)的協(xié)方差矩陣一致，因此，對(duì)于實(shí)際問題，這將導(dǎo)致耗費(fèi)大量計(jì)算資源。文獻(xiàn)[8]研究利用最小樣本數(shù)量進(jìn)行不確定度分析的方法，并提出了面向協(xié)方差矩陣的COST方法，但并未給出COST方法的理論基礎(chǔ)。為采用少樣本也能實(shí)現(xiàn)高精度多維隨機(jī)變量序列抽樣，本文首先證明COST方法的理論基礎(chǔ)，然后利用兩個(gè)變量的例題進(jìn)行驗(yàn)證分析，研究變量之間相關(guān)性的去除和施加的方法，最后進(jìn)行工程上的驗(yàn)證分析。

1 改進(jìn)的相關(guān)隨機(jī)變量模擬理論

設(shè)Xm×n=(X1,X2,…,Xn)m為需要產(chǎn)生的m組n維相關(guān)正態(tài)分布隨機(jī)變量序列，其均值u=(u1,u2,…,un)，協(xié)方差矩陣為A。產(chǎn)生n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量抽樣序列Ym×n=(Y1,Y2,…,Yn)，每個(gè)變量抽樣序列數(shù)為m。利用式(1)和獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量抽樣Yn即可產(chǎn)生隨機(jī)變量X：

X=ΣY+u

(1)

式中，矩陣Σ需要通過對(duì)協(xié)方差矩陣A分解得到。

V[X]=A=ΣΣT

(2)

由式(2)可知，僅將協(xié)方差矩陣A分解成矩陣Σ和ΣT相乘的形式，就可利用式(1)來(lái)產(chǎn)生多維相關(guān)正態(tài)分布隨機(jī)變量序列[5]。

上述理論是較嚴(yán)格的數(shù)學(xué)推導(dǎo)，但上述理論中有一個(gè)假設(shè)，即n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量抽樣序列Ym×n=(Y1,Y2,…,Yn)之間的協(xié)方差矩陣為單位矩陣I，在實(shí)際計(jì)算過程中，發(fā)現(xiàn)上述假設(shè)并不完全成立，特別是在抽樣樣本較少的情況下，n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量抽樣序列之間的相關(guān)性非常大。文獻(xiàn)[8]假設(shè)這個(gè)協(xié)方差矩陣為I*，且此協(xié)方差矩陣可分解為：

I*=MMT

(3)

由于n維隨機(jī)變量正態(tài)分布的相關(guān)性傳遞到了隨后的抽樣樣本上，導(dǎo)致實(shí)際產(chǎn)生的抽樣樣本的協(xié)方差矩陣并不完全是A，而是矩陣A與協(xié)方差矩陣I*的混合，由此導(dǎo)致了產(chǎn)生的n維相關(guān)隨機(jī)變量的協(xié)方差矩陣與真實(shí)的協(xié)方差矩陣并不能滿足完全一致。為能使得產(chǎn)生的n維相關(guān)隨機(jī)變量的協(xié)方差矩陣完全一致，首先需產(chǎn)生真正的n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布，假設(shè)真正的n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布為Y′，則Y和Y′的關(guān)系為：

Y′=M-1Y

(4)

式(2)就變?yōu)椋?/p>

V[X]=E[(X-E(X))(X-E(X))T]=

E[(ΣY′+u-u)(ΣY′+u-u)T]=

E[(ΣY′)(ΣY′)T]=E[ΣY′Y′TΣT]=

ΣV(Y′)ΣT=ΣΣT

(5)

改進(jìn)后的方法與之前的方法不同之處在于，首先需先求出n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)變量抽樣序列的協(xié)方差矩陣，然后利用式(3)和(4)將其轉(zhuǎn)換成真正的n維獨(dú)立標(biāo)準(zhǔn)正態(tài)分布隨機(jī)抽樣序列，在此基礎(chǔ)上把實(shí)際的相關(guān)性加上去，才能產(chǎn)生滿足協(xié)方差矩陣守恒的n維相關(guān)隨機(jī)變量序列。

本文提出的高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法的具體實(shí)現(xiàn)方法為：根據(jù)樣本總數(shù)N將概率區(qū)間等分為N個(gè)，使用RANDOM_NUMBER函數(shù)抽取隨機(jī)數(shù)ξ1，根據(jù)標(biāo)準(zhǔn)正態(tài)分布的分布函數(shù)與誤差函數(shù)的對(duì)應(yīng)關(guān)系計(jì)算ξ1對(duì)應(yīng)的概率，舍選抽樣直到保證N個(gè)等概率區(qū)間中每個(gè)區(qū)間均能抽取到有且僅有1個(gè)隨機(jī)數(shù)，則抽樣結(jié)束，即單變量拉丁超立方正態(tài)分布抽樣，重復(fù)上述抽樣可得到多序列抽樣樣本，然后利用式(4)進(jìn)行去除相關(guān)變量之間的相關(guān)性，即可得到高精度相關(guān)變量隨機(jī)數(shù)序列。

2 測(cè)試和驗(yàn)證

2.1 簡(jiǎn)單例題測(cè)試和驗(yàn)證

采用文獻(xiàn)[5]中的例題，要求產(chǎn)生隨機(jī)變量X1、X2、X3，其中X1、X2∈[-1,1]，X3∈[0,4]。3個(gè)隨機(jī)變量具有如下相關(guān)系數(shù)矩陣R：

(6)

利用Fortran2003語(yǔ)言分別開發(fā)了Cholesky和Jacobi分解法對(duì)矩陣R進(jìn)行計(jì)算，假設(shè)X1、X2、X3的標(biāo)準(zhǔn)偏差分別為0.1、0.2、0.3，均值為1、2、3。利用上述方法產(chǎn)生的20個(gè)樣本列于表1。采用高精度相關(guān)變量隨機(jī)數(shù)序列方法(改進(jìn)方法)與傳統(tǒng)的Jacobi方法和Cholesky方法的計(jì)算結(jié)果[5]亦列于表1，由表1可見，即使對(duì)于3個(gè)抽樣序列，傳統(tǒng)的Jacobi方法和Cholesky方法均需產(chǎn)生大量樣本(大于2萬(wàn))才能保證相關(guān)系數(shù)一致。而采用高精度相關(guān)變量隨機(jī)數(shù)序列方法，僅采用20個(gè)樣本，3個(gè)參數(shù)的相關(guān)系數(shù)即能保持與原相關(guān)系數(shù)一致(表2)。

表1 不同方法的相關(guān)系數(shù)對(duì)比Table 1 Comparison of correlation coefficient with different calculation methods

表2 利用改進(jìn)方法產(chǎn)生的20個(gè)樣本Table 2 Twenty samples using high-precision correlated variable sampling method

2.2 矩陣分解

為研究該方法對(duì)多個(gè)隨機(jī)數(shù)序列的適用性，本文對(duì)238U的輻射俘獲反應(yīng)道33群相關(guān)系數(shù)矩陣進(jìn)行了研究?；贓NDF/B-Ⅶ.1核評(píng)價(jià)數(shù)據(jù)庫(kù)，采用國(guó)際著名核截面處理程序NJOY2016，制作了238U的輻射俘獲反應(yīng)道33群相關(guān)系數(shù)矩陣，如圖1所示。由圖1可見，主對(duì)角線表示本能群與本能群的相關(guān)系數(shù)，所有能群的自相關(guān)系數(shù)均為1，并且不同能群之間具有較強(qiáng)的相關(guān)性。采用本文提出的高精度相關(guān)變量隨機(jī)數(shù)序列方法產(chǎn)生了33個(gè)抽樣序列，即n=33，每個(gè)抽樣序列樣本數(shù)為34，即m=34。研究發(fā)現(xiàn)，若m≤n，則式中的矩陣M無(wú)法求逆矩陣，因此，m的最小值為n+1。所以本文抽樣的每個(gè)抽樣序列的樣本數(shù)為34。圖2示出了利用抽樣得到的樣本序列計(jì)算得到的33個(gè)能群之間的相關(guān)系數(shù)矩陣，與圖1給出的原矩陣符合很好。圖3示出了抽樣樣本得到的相關(guān)系數(shù)矩陣(圖2)與原矩陣的相對(duì)誤差分布，最大相對(duì)誤差僅0.061%，絕大部分的相對(duì)誤差為0，由此說(shuō)明了高精度相關(guān)變量隨機(jī)數(shù)序列方法具有較強(qiáng)的處理多變量的能力，且具有高精度的特點(diǎn)。

圖1 238U 33群俘獲截面的相關(guān)系數(shù)矩陣Fig.1 Correlation coefficient matrix of 238U capture with 33 neutron groups

圖2 利用34個(gè)相關(guān)變量樣本序列計(jì)算得到的238U俘獲反應(yīng)道相關(guān)系數(shù)矩陣Fig.2 Calculated correlation coefficient matrix of 238U capture with 34 high-precision correlated variable pseudo-random number sequences

圖3 238U俘獲相關(guān)系數(shù)矩陣的相對(duì)誤差Fig.3 Relative error of correlation coefficient matrix of 238U capture

2.3 鉛基快堆LFR有效增殖因數(shù)不確定度

鉛基快堆作為世界上第四代先進(jìn)核反應(yīng)堆代表堆型之一，與傳統(tǒng)核反應(yīng)堆相比有著獨(dú)特的優(yōu)勢(shì)及特點(diǎn)，且有望最早實(shí)現(xiàn)商業(yè)化的第四代核反應(yīng)堆[9]。本文以歐共體核能局(NEA)組織的國(guó)際評(píng)價(jià)合作工作組(簡(jiǎn)稱NEAWPEC)提出的鉛基快堆LFR為研究對(duì)象[10]，研究主要核素不同反應(yīng)道對(duì)有效增殖因數(shù)keff的不確定度的貢獻(xiàn)，分析了高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法工程應(yīng)用的可行性。LFR的二維rz模型如圖4所示，該反應(yīng)堆主要由堆芯、反射層、屏蔽層、軸向屏蔽層、控制棒等組成。采用鉛(Pb)作冷卻劑、U-TRU-Zr合金作燃料，堆芯中次錒系核素(MA)的含量為2.0%，反應(yīng)堆滿功率為900 MWt，反應(yīng)堆循環(huán)壽期為310 d。

圖4 鉛基快堆LFR的rz模型Fig.4 rz model of lead cooled fast reactor LFR

選取了對(duì)LFR的keff具有重要貢獻(xiàn)的6個(gè)核素多個(gè)反應(yīng)道進(jìn)行研究，對(duì)比高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法與傳統(tǒng)抽樣方法計(jì)算得到keff的不確定度的相對(duì)誤差。239Pu和238U均為L(zhǎng)FR燃料中的重要核素，206Pb、207Pb和208Pb是本基準(zhǔn)例題LFR冷卻劑的成分，56Fe為結(jié)構(gòu)材料的重要組成部分，而燃料、冷卻劑和結(jié)構(gòu)材料幾乎構(gòu)成了LFR堆芯的全部材料，因此，這幾個(gè)核素具有很好的代表性。計(jì)算結(jié)果列于表3。表3中的參考解是采用直接擾動(dòng)方法得到每個(gè)核素每個(gè)反應(yīng)道的敏感性系數(shù)，然后再利用Sandwich公式計(jì)算得到。每個(gè)反應(yīng)道的抽樣樣本總數(shù)為50，傳統(tǒng)正態(tài)分布抽樣方法的抽樣樣本為100。由表3可見，雖然傳統(tǒng)正態(tài)分布抽樣樣本總數(shù)較高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法的樣本總數(shù)高1倍，但其相對(duì)誤差最大可達(dá)12.5%，而高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法的最大相對(duì)誤差僅1.7%，計(jì)算精度有明顯提高。

表3 不同計(jì)算方法計(jì)算各核素相對(duì)于keff的不確定度對(duì)比Table 3 Comparison of uncertainty of keff with respect to different isotopes by different sampling methods

3 結(jié)論

當(dāng)采用蒙特卡羅方法研究許多問題時(shí)，需要產(chǎn)生多維相關(guān)隨機(jī)數(shù)序列。但之前的研究表明：需特別大的抽樣樣本序列才能基本保證樣本的協(xié)方差矩陣與真實(shí)的協(xié)方差矩陣一致。由于采用大量抽樣樣本，導(dǎo)致很大的時(shí)間消耗。本文從理論上分析了以前的多維相關(guān)變量隨機(jī)數(shù)序列的協(xié)方差矩陣與真實(shí)的協(xié)方差矩陣有差別的原因是之前產(chǎn)生的多維相關(guān)變量隨機(jī)數(shù)序列，在考慮真實(shí)協(xié)方差矩陣前，相關(guān)變量隨機(jī)數(shù)序列之間并不能完全滿足不相關(guān)假設(shè)條件，通過對(duì)產(chǎn)生后的樣本的進(jìn)一步處理，即可得到，在考慮真實(shí)協(xié)方差矩陣之前，即能滿足各變量之間完全滿足不相關(guān)條件。然后，利用數(shù)值方法對(duì)上述理論進(jìn)行了驗(yàn)證，通過3個(gè)變量的例題可見，原來(lái)的抽樣方法需上萬(wàn)個(gè)抽樣樣本才能計(jì)算得到與原相關(guān)協(xié)方差矩陣一致的矩陣，而采用高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法，20個(gè)抽樣樣本即可完全重構(gòu)出原協(xié)方差矩陣。樣本數(shù)量減小了5個(gè)量級(jí)，而計(jì)算精度有了提高。在此基礎(chǔ)上，通過對(duì)33群的238U輻射俘獲反應(yīng)的相關(guān)系數(shù)矩陣進(jìn)行分析，得到了產(chǎn)生的樣本數(shù)一定要大于樣本變量數(shù)，利用高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法，產(chǎn)生的樣本可高度還原原相關(guān)系數(shù)矩陣，最大相對(duì)誤差小于0.061%。最后，利用鉛基快堆LFR對(duì)該方法的工程應(yīng)用性進(jìn)行了研究，研究結(jié)果表明，雖然傳統(tǒng)正態(tài)分布抽樣樣本總數(shù)較高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法的樣本總數(shù)高1倍，但其相對(duì)誤差仍較大，最大達(dá)12.5%，而高精度相關(guān)變量隨機(jī)數(shù)序列產(chǎn)生方法的最大相對(duì)誤差僅1.7%，計(jì)算精度有明顯提高。采用此方法時(shí)，抽樣序列樣本數(shù)要大于變量的個(gè)數(shù)，即m>n，否則，因?yàn)榫仃嘙無(wú)法求逆矩陣，而無(wú)法采用該方法，m的最小值為n+1。通過鉛基快堆LFR驗(yàn)證表明，該方法具有工程應(yīng)用前景。