基于貝葉斯MCMC方法的高斯煙羽模型不確定性分析

崔威杰 曹 博 陳義學(xué)

1（華北電力大學(xué)核科學(xué)與工程學(xué)院 北京102206）

2（華北電力大學(xué)非能動(dòng)核能安全技術(shù)北京市重點(diǎn)實(shí)驗(yàn)室 北京102206）

當(dāng)前核電廠采取了越來(lái)越好的安全設(shè)施，安全性已經(jīng)達(dá)到了很高的水平，其中三代核電站AP1000發(fā)生大規(guī)模放射性物質(zhì)泄漏事故的概率已被降低到1 × 10-8a-1［1］，但仍然不能就此不再考慮核電廠發(fā)生嚴(yán)重事故的可能性。一旦核電廠發(fā)生了放射性物質(zhì)泄漏事故，放射性物質(zhì)隨大氣的擴(kuò)散將會(huì)造成快速、大范圍的嚴(yán)重影響，因此必須及時(shí)評(píng)估放射性后果以制定合理的應(yīng)對(duì)策略。這種情況下，由于風(fēng)場(chǎng)、大氣湍流、溫度層結(jié)等的復(fù)雜性和難以準(zhǔn)確監(jiān)測(cè)的問(wèn)題［2］，造成了氣象觀測(cè)資料的不確定性。此外，大氣擴(kuò)散模型中還存在由于對(duì)實(shí)際問(wèn)題的簡(jiǎn)化而導(dǎo)致的不確定性。如果在事故后果評(píng)價(jià)時(shí)直接使用模型得到的確定性結(jié)果，則可能出現(xiàn)決策的失誤［3］。因此需要采取合適的方法對(duì)大氣擴(kuò)散模型進(jìn)行不確定性分析，降低決策失誤的風(fēng)險(xiǎn)。

高斯煙羽模型是當(dāng)前大氣擴(kuò)散模型中應(yīng)用最為廣泛的模型之一，其突出優(yōu)點(diǎn)有輸入?yún)?shù)少、所需計(jì)算資源少、比較符合觀測(cè)情況等。2018年發(fā)布的環(huán)境影響評(píng)價(jià)技術(shù)導(dǎo)則（HJ 2.2-2018）［4］中，推薦了兩種基于高斯煙羽模型的進(jìn)一步預(yù)測(cè)模型——AERMOD［5］和 ADMS［6］。除此之外 ，美國(guó)開(kāi)發(fā)的HotSpot［7］程 序 、歐 洲 開(kāi) 發(fā) 的 ATSTEP［8］模 型 及Polyphemus系統(tǒng)［9］，以及廣泛用于核電廠評(píng)審的ARCON96模型［10］和PAVAN程序［11］等都包含基于高斯煙羽模型的計(jì)算模塊。本文利用基于高斯煙羽模型及其修正條件開(kāi)發(fā)的大氣擴(kuò)散模擬程序進(jìn)行了參數(shù)敏感性分析和不確定性分析。

1 高斯煙羽模型及其敏感性分析

1.1 理論模型

本文使用高斯煙羽模型作為計(jì)算模型，以連續(xù)點(diǎn)源釋放引起的大氣放射性濃度作為計(jì)算對(duì)象。該模型假設(shè)放射性污染物釋放點(diǎn)在位于地面上的原點(diǎn)的正上方，下風(fēng)向?yàn)閤軸方向，y軸垂直于x軸，z軸由原點(diǎn)垂直向上。模型簡(jiǎn)化過(guò)程中采用了以下基本假設(shè)［12］：

1）污染物濃度在y軸、z軸方向上呈高斯分布（正態(tài)分布）；

2）風(fēng)速是均勻穩(wěn)定的，且大于1 m·s-1；

3）源強(qiáng)是均勻連續(xù)的點(diǎn)源；

4）湍流場(chǎng)是平穩(wěn)均勻的；

5）彌散過(guò)程中污染物質(zhì)量（或核素活度）是守恒的；

6）地面對(duì)煙羽具有完全反射作用。

基于以上假設(shè)，高斯煙羽模型的計(jì)算公式可以寫(xiě)成：

式中：Q是泄漏源強(qiáng)度，Bq·s-1；σy和σz分別為y方向和z方向上的高斯擴(kuò)散系數(shù)，m；u為有效釋放高度處的風(fēng)速，m·s-1；H為有效釋放高度，m。在上述方程的基礎(chǔ)上，考慮可導(dǎo)致煙羽耗減的干沉積、濕沉積以及放射性核素衰變等的影響后，可以得到更符合現(xiàn)實(shí)情況的高斯煙羽模型。

1.2 參數(shù)敏感性分析

在進(jìn)行模型參數(shù)不確定性分析之前，通常先進(jìn)行參數(shù)敏感性分析，以找出對(duì)模型結(jié)果影響較大的幾個(gè)參數(shù)，對(duì)這幾個(gè)參數(shù)進(jìn)行不確定性分析，可以在基本保證不確定性分析的準(zhǔn)確性和可信度時(shí)，大大減少計(jì)算量。

參數(shù)敏感性分析的主要方法有一次一個(gè)變量法、標(biāo)準(zhǔn)回歸系數(shù)法及Sobol敏感度指標(biāo)等方法。一次一個(gè)變量法是每次變動(dòng)一個(gè)參數(shù)，保持其余參數(shù)不變，統(tǒng)計(jì)運(yùn)算結(jié)果的變動(dòng)情況，進(jìn)而得出模型結(jié)果對(duì)該參數(shù)的敏感性；標(biāo)準(zhǔn)回歸系數(shù)法是用隨機(jī)抽樣方法產(chǎn)生輸入?yún)?shù)樣本，將其輸入到模型中得到模型結(jié)果，進(jìn)行最小二乘回歸分析；Sobol敏感度指標(biāo)按照參數(shù)對(duì)模型結(jié)果方差的貢獻(xiàn)來(lái)劃分敏感性等級(jí)［13］。

本文采用一次一個(gè)變量方法對(duì)5個(gè)重要輸入?yún)?shù)（釋放高度、10 m高度處風(fēng)速（以下稱(chēng)參考風(fēng)速）、煙羽釋放半徑、煙羽釋放速率及煙羽溫度）進(jìn)行敏感性分析。首先確定各參數(shù)初始值，然后以初始值的20%為步長(zhǎng)進(jìn)行依次計(jì)算，記錄每個(gè)參數(shù)變化過(guò)程中模型計(jì)算出的最大空氣濃度值，最后分析每個(gè)參數(shù)變化引起的最大空氣濃度標(biāo)準(zhǔn)差，標(biāo)準(zhǔn)差越大則表明模型對(duì)相應(yīng)參數(shù)的敏感性越大［14］。源項(xiàng)假設(shè)為單核素137Cs釋放，釋放率為3.70×1010Bq·s-1，大氣穩(wěn)定度設(shè)置為D。分析過(guò)程如表1所示。

由表1中的數(shù)據(jù)計(jì)算得到釋放高度、參考風(fēng)速、釋放半徑、釋放速率、煙羽溫度依次對(duì)應(yīng)的最大濃度標(biāo) 準(zhǔn) 差 依 次 為 1.09×105Bq·m-3、4.71×104Bq·m-3、2.17×104Bq·m-3、1.09×104Bq·m-3、8.56×103Bq·m-3，因此參數(shù)敏感性從大到小排序?yàn)椋横尫鸥叨龋瑓⒖硷L(fēng)速，煙羽釋放半徑，煙羽釋放速率，煙羽溫度。模型對(duì)釋放高度和參考風(fēng)速具有最大的敏感性，故而選取這兩個(gè)參數(shù)進(jìn)行不確定性分析。

表1 一次一個(gè)變量法計(jì)算過(guò)程Table 1 Calculation process of one-variable-at-a-time approach

2 貝葉斯MCMC方法

2.1 貝葉斯基本理論

貝葉斯理論最早由英國(guó)學(xué)者貝葉斯在18世紀(jì)中葉提出，在20世紀(jì)50年代后，貝葉斯理論得到了應(yīng)有的重視，并開(kāi)始在統(tǒng)計(jì)決策問(wèn)題中得到廣泛應(yīng)用［15］。貝葉斯理論的核心是貝葉斯公式，可以表示為：

式中：θ為模型輸入?yún)?shù)構(gòu)成的向量；p(θ)為這些參數(shù)的聯(lián)合先驗(yàn)概率密度；p(z|θ)為模型輸出值與實(shí)際觀測(cè)值z(mì)在輸入?yún)?shù)向量θ時(shí)的相似度，通常用似然函數(shù)表示。

2.2 MCMC方法

由于高斯煙羽模型經(jīng)修正后很難找到準(zhǔn)確的計(jì)算表達(dá)式，無(wú)法直接應(yīng)用貝葉斯公式求解參數(shù)不確定性，因此采用了馬爾科夫鏈蒙特卡羅方法（Markov Chain Monte Carlo，MCMC）來(lái)耦合貝葉斯方法和高斯煙羽模型［16］。常用的MCMC算法有Metropolis-Hastings算法、Gibbs抽樣和自適應(yīng)Metropolis算法等。下面介紹MCMC方法的自適應(yīng)Metropolis算法［17］：

1）按照先驗(yàn)分布隨機(jī)抽樣產(chǎn)生初始樣本向量θ0；

2）產(chǎn)生候選樣本向量：由期望為上個(gè)樣本向量θt-1及協(xié)方差矩陣為Ct經(jīng)隨機(jī)抽樣產(chǎn)生候選樣本向量x，其中：

式中：t0為初始化階段的樣本數(shù)；C0為初始協(xié)方差矩陣；COV(θ0，θ1，…，θt-1)表示所有已有的樣本向量的協(xié)方差；Id為維數(shù)是樣本向量維數(shù)d的單位矩陣；ε是一個(gè)很小的數(shù)，用于防止協(xié)方差降為0；sd是一個(gè)比例因子，用于保證合適的接受概率，此處sd=2.42/d［18］。

協(xié)方差矩陣可由式（4）計(jì)算：

3）產(chǎn)生樣本向量θt：以候選樣本向量作為模型的輸入?yún)?shù)，得到模型計(jì)算結(jié)果，然后以下面的接受概率確定是否接受候選樣本向量：

同時(shí)隨機(jī)抽樣取一個(gè)0～1之間的隨機(jī)數(shù)μ，若μ≤α，則θt=x，反之θt=θt-1。

4）重復(fù)步驟2）和步驟3），直到獲得足夠數(shù)量的參數(shù)樣本向量。

5）判斷樣本向量是否收斂到貝葉斯后驗(yàn)分布。

3 高斯煙羽模型不確定性分析

3.1 獲得觀測(cè)資料

由于核電廠放射性物質(zhì)泄漏事故的觀測(cè)資料難以獲得，因此本文使用大氣擴(kuò)散程序生成模擬數(shù)據(jù)作為觀測(cè)資料使用。將假設(shè)的理想氣象條件輸入到程序中，得到幾個(gè)觀測(cè)點(diǎn)處的大氣放射性濃度作為觀測(cè)資料。生成模擬觀測(cè)數(shù)據(jù)時(shí)，釋放高度設(shè)置為50 m，參考風(fēng)速設(shè)置為5.0 m·s-1。觀測(cè)點(diǎn)高度設(shè)置為1.5 m，共設(shè)置了5個(gè)觀測(cè)點(diǎn)，分別位于下風(fēng)向0.5 km、0.8 km、1.0 km、2.0 km、5.0 km處。

3.2 似然函數(shù)

在貝葉斯分析中，將模型輸入?yún)?shù)中的釋放高度和風(fēng)速作為隨機(jī)變量。似然函數(shù)的形式為：

式中：θ表示參數(shù)樣本向量；i=1，2，…，n，n為觀測(cè)點(diǎn)數(shù)；F是一個(gè)算子，表示將參數(shù)樣本輸入到模型中計(jì)算得到的指定觀測(cè)點(diǎn)的空氣放射性比濃度值；z*i表示第i個(gè)觀測(cè)點(diǎn)的空氣放射性比濃度觀測(cè)值。

3.3 MCMC方法應(yīng)用于高斯煙羽模型

3.3.1 確定先驗(yàn)分布

自適應(yīng)Metropolis算法相比于其他MCMC算法的最大優(yōu)勢(shì)在于不再需要嚴(yán)格確定參數(shù)的先驗(yàn)分布類(lèi)型，只需要確定輸入?yún)?shù)可能取值的范圍即可［19］，因此本文假設(shè)釋放高度和10 m高度處的風(fēng)速的先驗(yàn)分布為均勻分布，分布范圍如表2所示。

3.3.2 MCMC抽樣

先驗(yàn)分布確定后，初始協(xié)方差矩陣選為參數(shù)先驗(yàn)分布范圍的1/20，較小的初始協(xié)方差矩陣能使得MCMC樣本序列更快收斂。初始化階段樣本數(shù)為200，總的采樣樣本數(shù)為40 000。圖1和圖2給出了樣本均值和標(biāo)準(zhǔn)差隨迭代過(guò)程的變化曲線。

表2 輸入?yún)?shù)的先驗(yàn)分布Table 2 Prior distribution of input parameters

圖1 輸入?yún)?shù)樣本均值的變化曲線(a)釋放高度，(b)參考風(fēng)速Fig.1 Sample mean change curve of input parameters(a)Release height,(b)Wind speed

從上面的樣本序列均值和標(biāo)準(zhǔn)差變化曲線可以看出，在20 000次采樣之后，樣本序列的均值和標(biāo)準(zhǔn)差曲線的波動(dòng)幅度均較之前大大減小，即趨于穩(wěn)定。除上述判斷標(biāo)準(zhǔn)外，Gelman等［20］提出了基于多個(gè)采樣序列的收斂判斷方法：

式中：R為比例降低系數(shù)，它的值越接近1，說(shuō)明樣本序列越趨于收斂；n為每個(gè)樣本序列的采樣數(shù)；m為樣本序列數(shù)；B/n為樣本序列的均值之間的標(biāo)準(zhǔn)差；W為樣本標(biāo)準(zhǔn)差之間的均值。

圖2 輸入?yún)?shù)樣本標(biāo)準(zhǔn)差的變化曲線(a)釋放高度，(b)參考風(fēng)速Fig.2 Sample standard deviation change curve of input parameters(a)Release height,(b)Wind speed

為減少采樣偶然性的影響，并使用上述收斂判斷方法，重復(fù)進(jìn)行了3次MCMC采樣，每次獲取40 000個(gè)樣本。計(jì)算得出采樣過(guò)程中R的變化曲線，如圖3所示。

可以看出，在采樣初期，比例降低系數(shù)R變化劇烈，呈現(xiàn)波動(dòng)下降趨勢(shì)；在20 000次采樣以后，兩個(gè)輸入?yún)?shù)的R值都相當(dāng)接近1，這說(shuō)明此時(shí)MCMC序列已經(jīng)收斂。

3.4 樣本序列分析

從上一步的分析可知，樣本序列在20 000次采樣之后收斂，因此取每個(gè)樣本序列的后20 000個(gè)樣本，這樣就得到了60 000個(gè)收斂到貝葉斯后驗(yàn)分布的樣本及其對(duì)應(yīng)模擬結(jié)果。對(duì)其進(jìn)行分析得到每個(gè)觀測(cè)點(diǎn)處的大氣放射性濃度值的最小值、最大值、眾數(shù)以及標(biāo)準(zhǔn)差，如表3所示。

圖3 采樣過(guò)程中的比例降低系數(shù)Fig.3 The scale reduction score during sampling

表3中的眾數(shù)并不是傳統(tǒng)意義上出現(xiàn)次數(shù)最多的一個(gè)值，而是將樣本的取值范圍從小到大平分為50組后，取樣本點(diǎn)最多的一組的中間值作為眾數(shù)，這個(gè)值是模擬結(jié)果和觀測(cè)值相差最小的數(shù)值，可以認(rèn)為是觀測(cè)資料的最優(yōu)擬合。由于收斂到后驗(yàn)分布的樣本序列符合正態(tài)分布，因此模型結(jié)果也符合正態(tài)分布，統(tǒng)計(jì)模型結(jié)果的數(shù)值分布可以獲得置信區(qū)間。據(jù)此繪制出觀測(cè)資料的最優(yōu)擬合曲線以及模擬結(jié)果的95%置信區(qū)間如圖4所示。

圖4 觀測(cè)資料的最優(yōu)擬合及模擬結(jié)果的95%置信區(qū)間Fig.4 Optimal fitting and simulation results of observation data in 95%confidence interval

表3 模型結(jié)果序列分析Table 3 Analysis of model results series

4 結(jié)語(yǔ)

相比于傳統(tǒng)的不確定性分析方法，貝葉斯方法充分利用了已有的觀測(cè)資料和參數(shù)先驗(yàn)分布信息，其結(jié)果更加可靠，得到的置信區(qū)間的可信度也更高。馬爾科夫鏈蒙特卡羅方法可以方便地耦合大氣擴(kuò)散模型和貝葉斯方法，避免了貝葉斯公式中歸一化常數(shù)的求解，因此可以很好地應(yīng)用在難以尋求解析式的模型中。

根據(jù)MCMC方法樣本序列對(duì)應(yīng)的模型模擬結(jié)果的眾數(shù)和標(biāo)準(zhǔn)差數(shù)據(jù)分析，可以得到觀測(cè)值的最優(yōu)擬合和置信區(qū)間，從而為事故后應(yīng)急響應(yīng)提供更加可靠的數(shù)據(jù)參考，有效降低決策失誤的風(fēng)險(xiǎn)。