基于MCMC方法的子通道程序模型參數(shù)不確定性量化分析研究

何鑫 宋美琪 劉曉晶,2

1（上海交通大學 智慧能源創(chuàng)新學院 上海 200240）

2（上海交通大學 核科學與工程學院 上海 200240）

相較于保守性安全分析，最佳估算結合不確定性分析能更加真實準確地模擬反應堆運行狀況，更適合核電廠的安全裕度分析與性能評價［1］。1989年，美國核管理委員會發(fā)布了RGI.157管理導則，允許在應急堆芯冷卻系統(tǒng)（Emergency Core Cooling System，ECCS）分析中使用最佳估算程序，但要對計算結果的不確定性進行評估［2］。此后，多種與現(xiàn)實性程序相關的不確定性分析方法被開發(fā)出來，如著名的程序縮放模擬、應用及不確定性分析（Code Scaling，Applicability and Uncertainty，CSAU）方法［3］、西屋公司采用的ASTRUM（Advanced Statistical Treatment of Uncertainty Method）方法［4］、德國開發(fā)的GRS（Gesellschaft fur Anlagen-und Reaktorsicherheit）方法［5］、意大利開發(fā)的UMAE（Uncertainty Method based on Accuracy Extrapolation）方法［6］等，主要可以分為基于輸入的不確定性正向傳遞和基于輸出的不確定性外推，前者的發(fā)展更為成熟，應用更加廣泛［7］。這些方法大多要提前給定不確定性來源的分布，但受制于數(shù)值模型近似和認知水平局限，大多數(shù)熱工水力程序并沒有提供足夠的與內部模型相關的輸入不確定性信息，通常由專家評估給出。

近年來，一些代替專家保守評價的分析方法陸續(xù)被提出，如利用快速傅里葉變化來評估輸入?yún)?shù)變化范圍的FFTBM（Fast Fourier Transformation Based Methods）方法［8］、通過不斷調整參數(shù)邊界來評估程序計算值對實驗數(shù)據(jù)的覆蓋率得到參數(shù)的偽累積分布函數(shù)的DIPE（Determination of Input Parameters Empirical properties）方法［9］、基于線性假設通過EM（Expectation Maximization）算法評估模型參數(shù)不確定性的CIRCé（Calcul des Incertitudes Relatives aux Corrélations élémentaires）方法［10］，以及通過確定性和概率性數(shù)據(jù)同化方法求解線性和非線性問題的MCDA（Model Calibration through Data Assimilation method）方法［11］等。Wilks 等［12］利用Metropolis 算法實現(xiàn)馬爾科夫鏈蒙特卡羅（Markov Chain Monte Carlo，MCMC）抽樣，對收斂到平穩(wěn)目標分布的馬爾科夫鏈進行抽樣統(tǒng)計，從而獲得反問題的解。隨機論對應著大量的正向計算，因此建立高精度的代理模型來替代復雜的正向程序有著重要意義。李冬等［13］對MCDA 方法的非線性部分進行改進，通過建立高精度的RBF（Radial Basis Function）神經網(wǎng)絡代理模型來簡化正向計算，并將其應用于再淹沒現(xiàn)象中物理模型參數(shù)不確定性的評估。Wu 等［14］采用基于高斯過程（Gaussian Process，GP）代理模型的改進模塊化貝葉斯方法對TRACE程序的相關物理模型參數(shù)進行不確定性分析。

子通道程序是重要的堆芯熱工水力安全分析工具，但由于對熱工水力模型的認識仍存在欠缺，在相關的設計和安全分析中會采用一些保守估計，導致程序計算結果存在一定不確定性。而現(xiàn)階段的不確定性分析大多針對系統(tǒng)程序，有關子通道程序模型參數(shù)的不確定性量化分析相對較少。本研究采用隨機論方法，利用MCMC抽樣構造以模型參數(shù)不確定性后驗分布為平穩(wěn)分布的馬爾科夫鏈進行抽樣統(tǒng)計，并通過相對熵最小化自適應加密訓練數(shù)據(jù)從而構建高精度BPNN（Back Propagation Neural Network）代理模型代替復雜的系統(tǒng)程序，從而對影響空泡份額的模型參數(shù)進行了不確定性量化分析。

1 模型參數(shù)的不確定性量化分析方法

不確定性分析是指識別和量化影響輸出的所有不確定性因素，特別是不能通過實驗直接測量的關鍵模型參數(shù)，并基于不確定性的正向傳遞得到輸出結果總的不確定性分布［7］。不確定性的來源［15］有程序模型偏差，數(shù)值計算近似，邊界條件和初始條件不充分等，其中程序模型對物理現(xiàn)象的近似是不確定性的根本來源。

1.1 模型參數(shù)的不確定性量化

不確定性分析需要確定模型參數(shù)基準值的無量綱常數(shù)k（k=(kj)j=1...p，p指模型參數(shù)個數(shù)），可用kj=exj表示［16］。其中，x=(xj)j=1...p指模型參數(shù)的不確定性，y=(yj)j=1..q（q指實驗點個數(shù)），表示程序響應的不確定性，不確定性在正向程序中的傳遞關系可簡化為y=h(x)。假設參數(shù)之間相互獨立，可以用正態(tài)分布來描述未知模型參數(shù)的不確定，即x～N(μx，Cx)。實驗測量值yexp與程序預測值ycode為真實值yreal的求解提供了參考，兩者差值ym可用式（1）表示，其中：e指實驗誤差，且e～N(0，σexp)。

根據(jù)貝葉斯原理，模型參數(shù)的后驗概率ppost(x|ym)可表示為式（2），其中pprior(x)是模型參數(shù)不確定性的先驗概率密度，p(ym)指可觀測誤差的概率，可看作常數(shù)。L(ym|x)指模型參數(shù)不確定性x關于可觀測誤差ym的似然函數(shù)。

本研究通過Metropolis算法實現(xiàn)MCMC抽樣求解式（2）。首先，需要構建一個合理的提議函數(shù)，并根據(jù)接受概率α判斷是否保留第t步的采樣點x*，x*在x*～N(xt-1，Ct-1)中抽樣產生，直至通過循環(huán)產生的點構成的馬爾科夫鏈可以穩(wěn)定收斂到后驗分布。接受概率表示為式（3）。

然而，要想得到穩(wěn)定的馬爾科夫鏈，抽樣次數(shù)至少達到104數(shù)量級，在求解公式（2）時需要調用相應次數(shù)的最佳估算程序，這意味著巨大的計算量。為了提高計算效率，可通過BPNN 建立代理模型代替式（2）中的

1.2 BPNN 代理模型構建及訓練數(shù)據(jù)自適應加密

BPNN 對非線性函數(shù)的擬合性較好，具有收斂速度快，魯棒性好，適應性強的優(yōu)點。BPNN是利用誤差反向傳播算法進行訓練的多層前饋神經網(wǎng)絡，它將網(wǎng)絡輸出的誤差平方作為目標函數(shù)，利用梯度下降法求解目標函數(shù)的最小值實現(xiàn)權值的更新。其中輸入層的樣本點{x1，x2，...，xn}是在上一步迭代樣本空間中通過隨機抽樣產生，隱藏層采用sigmoid函數(shù)，輸出層通過運行外部正向程序得到x對應的構成。

為提高代理模型的精度，訓練樣本點要盡可能在接近后驗概率的空間上取到，可利用相對熵最小化在高后驗概率空間自適應加密樣本點［17］。提議空間ppro(x)與后驗概率空間ppost(x|ym)的差異表示為式（4），其展開式的前一項與ppro(x)無關。因此最小化DKL(ppost||ppro)等同于最大化后一項，記為D(θpro)，表示為式（5）?？赏ㄟ^來不斷更新提議分布的均值與方差，直至滿足收斂要求。在建立滿足精度要求的代理模型后，給定隨機抽樣產生的輸入樣本點x，就可以利用BPNN代理模型得到對應的輸出。

1.3 不確定性分析方法驗證

本文設置了一個簡單算例驗證上述方法的準確性，其函數(shù)表達式如式（6）所示，其中，x1與x2相互獨立，可知其最大概率點位于點（3，4），方差分別是0.4和0.3。

圖1 不確定性分析流程圖Fig.1 Flowchart of uncertainty analysis

表1 自適應加密過程中的提議參數(shù)Table 1 Proposed parameters for adaptive encryption

2 子通道程序不確定性分析

2.1 子通道程序對PSBT 空泡份額實驗不確定性量化分析

壓水堆子通道及管束試驗（PWR Sub-channel and Bundle Tests，PSBT）是OECD/NRC 發(fā)布的一套基準，主要用于驗證子通道和管束中空泡份額及偏離泡核沸騰（Departure from Nucleate Boiling，DNB），實驗數(shù)據(jù)包括單相流管束出口的溫度分布、穩(wěn)態(tài)和瞬態(tài)的空泡份額以及臨界功率［19］?？张莸拇嬖跁绊懛磻讯研景羰膫鳠嵝?，所以能否高效準確地預測空泡份額具有重要意義。

本研究通過子通道程序COBRA-IV 模擬PSBT基準一階段穩(wěn)態(tài)空泡份額分布，從而對子通道程序的計算能力和不確定性進行分析。選取PSBT穩(wěn)態(tài)工況的30組空泡份額實驗數(shù)據(jù)，其壓力變化范圍是7.36～14.71 MPa，質量流量變化范圍是1 380.6～2 288.9 kg·（m2·s）-1，加熱功率變化范圍是1 920～3 536 kW，入口溫度變化范圍是173.5～287.8 ℃。選取其中的25組，利用開發(fā)的不確定性分析程序反向量化對空泡份額有重要影響的模型參數(shù)的不確定性，在確定模型參數(shù)的不確定性后，基于不確定性正向傳遞的方法，利用另外5 組實驗數(shù)據(jù)對不確定性結果進行評估。

在選用滑移模型作為空泡份額模型時，滑速比和湍流交混系數(shù)是影響空泡份額的主要模型參數(shù)［20］，其他影響空泡份額的邊界條件參數(shù)不確定性可認為服從均值為0 的正態(tài)分布，標準差可通過PSBT 實驗測量精度確定，并按照3σ原則進行處理［21］，結果如表2 所示。由于多波束系統(tǒng)測量弦向平均孔隙率的穩(wěn)態(tài)精度為±4%，可確定實驗誤差σexp為4%［21］。

表2 邊界條件參數(shù)的不確定性分布Table 2 Uncertainty distribution of boundary condition parameters

本研究采用滑移模型，其滑速比可以進行設置?；俦群屯牧鹘换煜禂?shù)根據(jù)如下標準進行選?。?2］：根據(jù)最小熵增原理可計算14 MPa 下的滑移比為1.86，但通常最佳估算程序的計算值要大于實驗測量值；均相模型是最典型的單向流模型，認為氣相和液相之間不存在相對滑移，即滑速比為1，因此，最終將滑速比的基準值設為1.2。Castellana實驗與本研究選擇的工況相似，因此取湍流交混系數(shù)的基準值為0.008。

將滑速比和湍流交混系數(shù)設置為選定的基準值進行計算，子通道程序的空泡份額預測值和實驗測量值的對比如圖2 所示，數(shù)據(jù)點越靠近直線表明程序預測結果與實驗結果越接近，兩條虛線代表實驗誤差范圍?？梢园l(fā)現(xiàn)整體的預測值要高于實驗測量值，尤其是在低空泡份額有很多數(shù)據(jù)點落在了實驗誤差范圍外，這說明所選用的模型不確定性較大，需要對其進行量化分析。

圖2 空泡份額實驗值與計算值對比圖Fig.2 Comparison between the experimental and calculated void fraction values

2.2 模型參數(shù)不確定性分析結果

根據(jù)圖1 的不確定性分析流程圖，首先需要確定初次迭代的先驗分布，這里兩個模型參數(shù)服從相同的先驗分布，即sprior～N(0，0.1)，βprior～N(0，0.1)，其中：s為滑速比，β為湍流交混系數(shù)。每次迭代的抽樣個數(shù)為10，選取其中25組實驗工況用于模型參數(shù)的不確定性分析。將每組工況距離加熱底端2 669 mm 和3 177 mm 的中心子通道空泡份額平均值 作 為 實 驗 測 量 值 ，建 立g(x) =的代理模型。

利用相對熵最小化原理在自適應加密過程中得到的模型參數(shù)的提議分布如表3所示，可以看到，在三次迭代后滿足了收斂準則的前一項條件。由于g(x)表示50個樣本點的累計誤差，以及在計算中還存在著其他的隨機誤差，模型誤差等因素，所以會導致程序計算結果與實驗值不會完全重合。采用小批量隨機梯度下降法來對BPNN 進行訓練，在建立BPNN代理模型過程中每次迭代的網(wǎng)絡誤差下降曲線分別如圖3（a～c）所示。本研究嘗試了多種不同的初始先驗分布以及每次迭代時采用不同的抽樣點數(shù)，最終得到的提議參數(shù)與表3類似。

圖3 建立BPNN代理模型時的誤差下降曲線Fig.3 Error decline curve when establishing the BPNN surrogate model

表3 自適應加密迭代過程中得到的提議參數(shù)Table 3 Parameters proposed by adaptive algorithm iterations

在得到滿足要求的代理模型后，進行MCMC抽樣，設置鏈長為20 000。為了避免馬爾科夫鏈前期發(fā)散，對第5 000～20 000個點進行了統(tǒng)計分析，模型參數(shù)不確定性的結果如表4 所示。利用自適應Metropolis 算法得到的滑速比和湍流交混系數(shù)的隨機樣本游走點及頻次分布直方圖分別如圖4 和圖5所示。

圖4 自適應Metropolis 算法產生的樣本路徑圖 (a) s的樣本隨機游走，(b) β的樣本隨機游走Fig.4 Sample path of the adaptive Metropolis algorithm (a) Random-walk samples from s, (b) Random-walk samples from β

圖5 20 000次迭代樣本的頻次直方圖 (a) s的頻次直方圖，(b) β的頻次直方圖Fig.5 Frequency histogram of 20 000 iteration samples (a) Frequency histogram of s, (b) Frequency histogram of β

表4 模型參數(shù)不確定性分布的均值與標準差Table 4 Uncertainty distribution of the mean and standard variance for the model parameters

通過上述方法得到模型參數(shù)的不確定性后，同時考慮模型參數(shù)和邊界條件的不確定性，利用另外5 組實驗數(shù)據(jù)，基于不確定性的正向傳遞方法得到計算結果95%的置信區(qū)間對實驗值的包絡性情況。根據(jù)Wilks 公式，需要最少93 次抽樣才能得到雙邊容忍限和置信度均為95%的情況［23］，可通過蒙特卡羅方法進行不確定性的正向傳遞，運行93次最佳估算程序，將計算結果的最值作為置信區(qū)間的上下限。圖6 為雙95%置信區(qū)間的上下限包絡性檢驗，可以發(fā)現(xiàn)，不確定帶對實驗值的包絡性較好。利用模型參數(shù)不確定性統(tǒng)計均值對參數(shù)基準值進行修正，圖7為模型修正值和基準值的計算結果和實驗數(shù)據(jù)的對比，可以看出，修正后的模型預測值較原基準值預測結果更加接近實驗測量值，其中模型修正值計算結果與實驗值的平均絕對誤差（Mean Absolute Error，MAE）和均方根誤差（Root Mean Square Erroe，RMSE）分別為3.356%和4.079%，而基準值計算結果與實驗值的MAE 與RMSE 分別為12.028%和12.695%。

圖6 95%置信區(qū)間的包絡性檢驗 (a) 2 669 mm測量點的空泡份額，(b) 3 177 mm測量點的空泡份額Fig.6 Envelop test for 95% confidence interval (a) Void fraction at 2 669 mm, (b) Void fraction at 3 177 mm

圖7 模型參數(shù)均值修正后結果 (a) 2 669 mm測量點的空泡份額，(b) 3 177 mm測量點的空泡份額Fig.7 Calibration resultsof the two modified model parameters (a) Void fraction at 2 669 mm, (b) Void fraction at 3 177 mm

3 結語

本研究利用開發(fā)的Python 不確定性分析工具同時對多個模型參數(shù)不確定性進行分析，得到了滑速比和湍流交混系數(shù)的不確定性分布。在所選的計算模型和實驗工況范圍內，滑速比的不確定性較湍流交混系數(shù)更大。得到輸入?yún)?shù)不確定性后，基于不確定性正向傳遞方法，得到了95%置信區(qū)間的空泡份額不確定帶，并利用模型參數(shù)不確定性均值對基準值進行修正，結果表明：

1）原模型基準值預測結果要普遍高于實際測量值，尤其在低空泡份額情況，很多點落在了精度范圍外?；诓淮_定性正向傳遞方法，包絡性檢驗得到的95%置信區(qū)間不確定帶對實驗值包絡性較好。

2）利用模型參數(shù)不確定性均值對基準值進行修正，修正后的模型預測結果較基準值預測結果更加接近實驗值。因此本研究建立的不確定性量化分析方法能較好適用于子通道程序模型參數(shù)不確定性分析，這為后續(xù)的子通道程序不確定性研究提供了參考。

作者貢獻聲明何鑫負責研究方案設計，論文撰寫，實施研究，代碼編寫，采集數(shù)據(jù)，分析／解釋數(shù)據(jù)；宋美琪負責研究方案設計，對文章的知識性內容作批評性審閱，論文指導與修改，代碼調試，實驗結果數(shù)據(jù)分析；劉曉晶負責研究方案設計，獲取研究經費，技術支持，論文指導與修改，實驗結果數(shù)據(jù)分析，支持性貢獻等。

1D'Auria F, Camargo C, Mazzantini O. Best Estimate Plus Uncertainty (BEPU) approach in licensing current nuclear reactors[J]. Nuclear Engineering and Design, 2012, 248:317-328. DOI: 10.1016/j.nucengdes.2012.04.002.

3Skorek T, de Crécy A, Kovtonyuk A,et al. Quantification of the uncertainty of the physical models in the system thermal-hydraulic codes - PREMIUM benchmark[J].Nuclear Engineering and Design, 2019, 354: 110199.DOI: 10.1016/j.nucengdes.2019.110199.