崔威杰 曹 博 陳義學

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

2（華北電力大學非能動核能安全技術(shù)北京市重點實驗室 北京102206）

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

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

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

1.1 理論模型

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

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

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

3）源強是均勻連續(xù)的點源；

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

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

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

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

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

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

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

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

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

由表1中的數(shù)據(jù)計算得到釋放高度、參考風速、釋放半徑、釋放速率、煙羽溫度依次對應(yīng)的最大濃度標 準 差 依 次 為 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ù)敏感性從大到小排序為：釋放高度，參考風速，煙羽釋放半徑，煙羽釋放速率，煙羽溫度。模型對釋放高度和參考風速具有最大的敏感性，故而選取這兩個參數(shù)進行不確定性分析。

表1 一次一個變量法計算過程Table 1 Calculation process of one-variable-at-a-time approach

2 貝葉斯MCMC方法

2.1 貝葉斯基本理論

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

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

2.2 MCMC方法

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

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

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

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

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

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

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

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

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

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

3.1 獲得觀測資料

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

3.2 似然函數(shù)

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

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

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

3.3.1 確定先驗分布

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

3.3.2 MCMC抽樣

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

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

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

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

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

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

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

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

3.4 樣本序列分析

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

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

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

圖4 觀測資料的最優(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é)語

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

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