關(guān)鍵詞：放射性泄漏源參數(shù)反演;環(huán)境監(jiān)測數(shù)據(jù);大氣擴(kuò)散模型;迭代優(yōu)化;貝葉斯推斷

0 引言

近些年來，未知來源的放射性泄漏事件屢次發(fā)生，比如在2016 年1 月6 日，朝鮮開展了第四次地下核試驗(yàn)，一些惰性氣體監(jiān)測站臺(tái)探測到異常的133Xe 濃度，研究表明釋放位置為豐溪里核試驗(yàn)場[1] 。在2017 年10 月，歐洲大部分國家都報(bào)道了大氣中106Ru 濃度的升高，大部分研究表明最可能的釋放來源位于俄羅斯境內(nèi)的Mayak 核設(shè)施[2] 。為了評估放射性物質(zhì)在大氣中傳輸?shù)暮蠊?，核?yīng)急響應(yīng)系統(tǒng)基于氣象數(shù)據(jù)和大氣擴(kuò)散模型進(jìn)行模擬預(yù)報(bào)。但是當(dāng)泄漏源位置未知時(shí)，如何利用有限的環(huán)境和氣象信息，從空間和時(shí)間上量化泄漏源是一項(xiàng)極具挑戰(zhàn)的工作。這種利用環(huán)境監(jiān)測數(shù)據(jù)反推放射性泄漏源的位置和釋放率的 過程被稱為放射性泄漏源參數(shù)反演，是當(dāng)前核應(yīng)急研究的熱點(diǎn)問題。本文將討論分析放射性泄漏源參數(shù)反演系統(tǒng)的組成，論述反演方法的原理，綜述反演方法的研究成果，并在綜合比較方法性能的基礎(chǔ)上論述未來的發(fā)展方向。

1 放射性泄漏源參數(shù)反演系統(tǒng)介紹

泄漏源參數(shù)反演是源位置和釋放率的多解組合優(yōu)化過程，是大氣擴(kuò)散的逆向過程。如圖1 所示，放射性泄漏源參數(shù)反演系統(tǒng)包含四大要素：監(jiān)測數(shù)據(jù)、先驗(yàn)信息、大氣擴(kuò)散模型和反演方法。

1. 1 監(jiān)測數(shù)據(jù)

監(jiān)測數(shù)據(jù)常用空氣活度濃度數(shù)據(jù)，因?yàn)榭梢垣@取單一核素的數(shù)據(jù)，且具有較低的探測下限?！度娼购嗽囼?yàn)條約》（CTBT）要求建立的全球監(jiān)測系統(tǒng)包括80 個(gè)放射性核素監(jiān)測站臺(tái)[3] ，這些監(jiān)測站臺(tái)配備了氣溶膠采樣器和高精度的鍺探測器，其中40 個(gè)還額外配備了氙探測器，這些設(shè)備可以對空氣中的放射性核素進(jìn)行采樣和分析。對于某些放射性泄漏事件，由于信息封鎖，無法獲得靠近泄漏源的監(jiān)測數(shù)據(jù)，而且隨著核素在大氣中擴(kuò)散的傳輸，在某些時(shí)刻只有少數(shù)幾個(gè)站臺(tái)可以收集到核素樣品，因此監(jiān)測數(shù)據(jù)在時(shí)空上是稀疏的。

1. 2 先驗(yàn)信息

先驗(yàn)信息主要包含氣象場和先驗(yàn)估計(jì)，氣象場驅(qū)動(dòng)大氣擴(kuò)散模型運(yùn)行，主要的氣象來源包括歐洲中期天氣預(yù)報(bào)中心和美國國家環(huán)境預(yù)測中心，可分為再分析數(shù)據(jù)和預(yù)報(bào)數(shù)據(jù)。再分析數(shù)據(jù)是指利用過去的氣象觀測和模式預(yù)報(bào)結(jié)果，通過一套統(tǒng)一的氣象學(xué)模型和數(shù)據(jù)同化方法，對過去的氣象狀態(tài)進(jìn)行重新分析和再構(gòu)建，更接近真實(shí)情況;而預(yù)報(bào)數(shù)據(jù)則利用當(dāng)前的氣象觀測和模式預(yù)報(bào)結(jié)果，對未來某個(gè)時(shí)間段內(nèi)的氣象情況進(jìn)行推測，是不確定的，但是可用于預(yù)測未來的濃度分布。先驗(yàn)估計(jì)是根據(jù)現(xiàn)有資料對源參數(shù)初步的估計(jì)，需要一定的專家知識(shí)。對于核事故、核武器試驗(yàn)以及核材料泄漏等事件通常伴隨著大量放射性裂變產(chǎn)物的釋放，通過測量和分析空氣樣品中放射性裂變產(chǎn)物的同位素比率（如135Xe/ 133mXe，87Kr/ 88Kr，133I/ 131I 和137Cs/ 134Cs 等），可以大致推斷出放射性物質(zhì)首次被釋放的時(shí)刻，即核裂變反應(yīng)開始的時(shí)間，又稱之為核事件的零時(shí)刻，從而提供源釋放時(shí)間的先驗(yàn)信息[4] 。此外，可以利用非零監(jiān)測和零值監(jiān)測的反向軌跡確定釋放源大致區(qū)域，提供源位置的先驗(yàn)信息[5] 。非零監(jiān)測的反向軌跡途經(jīng)的位置更可能是實(shí)際泄漏源，通過疊加多個(gè)非零監(jiān)測的反向軌跡，可以初步計(jì)算源位置概率的空間分布;而零值監(jiān)測途經(jīng)的位置則可能性較小，可用于進(jìn)一步縮小源位置的可能范圍。然而，大多數(shù)反演方法通常僅假定源參數(shù)的大致范圍，并未設(shè)計(jì)專門的先驗(yàn)估計(jì)方法。

1. 3 大氣擴(kuò)散模型

放射性物質(zhì)泄漏后在大氣中主要經(jīng)過4 個(gè)過程：輸運(yùn)、擴(kuò)散、遷移和衰變，其濃度變化可由下式描述[6] ：

大氣擴(kuò)散模型是基于等式（1）進(jìn)行設(shè)計(jì)的集成化軟件，可分為兩種模式：拉格朗日模式和歐拉模式。拉格朗日模式假設(shè)核素以離散質(zhì)點(diǎn)形式存在，并隨大氣運(yùn)動(dòng)而移動(dòng)，通過求解質(zhì)點(diǎn)的軌跡來分析核素?cái)U(kuò)散的路徑和濃度。而歐拉模式則在固定的網(wǎng)格中計(jì)算核素濃度的變化，以直接求解公式（1）來模擬擴(kuò)散過程[7] 。在拉格朗日模式的基礎(chǔ)上，忽略復(fù)雜的湍流和不均勻擴(kuò)散現(xiàn)象，假設(shè)污染物擴(kuò)散過程中遵循高斯分布，可推導(dǎo)出高斯模式。表1 比較了這三種模式的優(yōu)缺點(diǎn)及其適用場景，并列舉了各模式的代表性模型。為結(jié)合拉格朗日模式和歐拉模式的優(yōu)勢，特別是在需要模擬局部與大尺度相結(jié)合的污染物傳輸時(shí)，能夠更好地平衡計(jì)算精度和效率，亦有采用拉格朗日-歐拉混合模式，如DREAM[14] 和HYPACT[15] 等。

如圖2 所示，大氣擴(kuò)散模型包括正向模式和反向模式（又稱伴隨模式）。正向模式是從放射性泄漏源（S）出發(fā)，模擬源向受體（R）的擴(kuò)散過程;反向模式則從受體出發(fā)，沿時(shí)間和空間反向追溯泄漏源。兩種模式均可計(jì)算出源-受體靈敏度（sourcereceptorsensitivity， SRS），并且現(xiàn)有研究表明，正向和反向模式的SRS 結(jié)果差異可忽略不計(jì)[16] 。顯然，當(dāng)潛在源的數(shù)量大于受體數(shù)量時(shí)，反向模式在計(jì)算上更加高效。這也符合反演計(jì)算的典型情況，即監(jiān)測點(diǎn)分布稀疏且源信息極其不明確。

由于大氣擴(kuò)散過程的復(fù)雜非線性以及氣象數(shù)據(jù)中的誤差，模擬結(jié)果與實(shí)際觀測數(shù)據(jù)之間不可避免地存在差異，這也導(dǎo)致了放射性泄漏源參數(shù)反演結(jié)果的誤差。為提高反演精度或更好地反映反演過程中存在的不確定性，一種有效的方法是采用集合模擬數(shù)據(jù)。具體而言，通過引入來自不同來源的氣象數(shù)據(jù)，或?qū)我粴庀髼l件進(jìn)行擾動(dòng)，構(gòu)造多組模擬數(shù)據(jù)，并基于這些數(shù)據(jù)進(jìn)行反演，獲得多組反演結(jié)果。隨后，對這些反演結(jié)果進(jìn)行統(tǒng)計(jì)分析，量化其不確定性，并采用中位數(shù)或均值等統(tǒng)計(jì)特征作為最終反演結(jié)果，從而提高反演的穩(wěn)健性[17] 。

1. 4 反演方法

放射性泄漏源可視為單一固定點(diǎn)源，假設(shè)泄漏源參數(shù)向量s = （r，q） 由源位置r = （lon，lat） 和釋放率q∈? T 成， T 表示時(shí)間步個(gè)數(shù)， M 個(gè)環(huán)境監(jiān)測數(shù)據(jù)組成向量μ ∈ ? M 。由于大氣擴(kuò)散模型誤差和監(jiān)測誤差的存在，即使采用真實(shí)源參數(shù)得到的模擬結(jié)果也無法與監(jiān)測值完全吻合，因此反演方法的本質(zhì)是最小化模擬結(jié)果與監(jiān)測值之間的差異[18] ：

式中， H ∈ ? M×T 是由SRS 組成的線性算子，每個(gè)元素表示源與某一受體之間的線性關(guān)系，但H 與r成非線性關(guān)系。由于監(jiān)測數(shù)據(jù)稀疏，等式（2）屬于非線性欠定問題，無法直接求解，因此反演方法一般基于帶有約束項(xiàng)的代價(jià)函數(shù)， 常用的形式如下[19] ：

該代價(jià)函數(shù)由兩項(xiàng)構(gòu)成，第一項(xiàng)描述模擬和監(jiān)測之間的差異，第二項(xiàng)為Tikhonov 正則項(xiàng)，用于給釋放率添加約束， λ 表示正則化參數(shù)，可由廣義交叉驗(yàn)證或者L-curve 等方法選取[20] 。

2 反演方法分類

根據(jù)放射性泄漏源參數(shù)反演的求解框架，反演方法可分為迭代優(yōu)化方法和貝葉斯推斷方法。假設(shè)釋放率q 的時(shí)間分辨率為ΔT ，其在不同的釋放行為假設(shè)下可簡化為不同的參數(shù)（釋放率曲線如圖3 所示）：

①瞬時(shí)釋放：釋放時(shí)間trel 和釋放量Q ;

②矩形釋放：釋放起始時(shí)間tstart ，釋放結(jié)束時(shí)間tend 和釋放率q ;

③ 變速釋放： 各個(gè)時(shí)間步內(nèi)的釋放率

2. 1 迭代優(yōu)化方法

迭代優(yōu)化方法首先遍歷所有可能源位置rn ，對每個(gè)源位置構(gòu)造H r（ n ） ，輸入監(jiān)測數(shù)據(jù)和H r（ n ） 到代價(jià)函數(shù)中，求得 qn = argminJ r（ n，q） ，找到使得J 最小的（ rn ，qn ） 即為反演結(jié)果。代價(jià)函數(shù)的形式多樣，最早有學(xué)者基于FLEXPART 模型和等式（3）對ETEX 實(shí)驗(yàn)的釋放源進(jìn)行反演，該方法非常適用于CTBT 的核素監(jiān)測系統(tǒng)[21] ;也有學(xué)者使用簡單的平方差代價(jià)函數(shù)μ - H rn （ ） q 22用于評估第四次朝核試驗(yàn)的泄漏源位置和釋放率[1] ;美國國家海洋和大氣管理局下屬的空氣資源實(shí)驗(yàn)室進(jìn)一步對模型不確定性進(jìn)行了建模，成功反演出CAPTEX 六次實(shí)驗(yàn)的源位置和釋放率結(jié)果[22] ;捷克科學(xué)院信息理論與自動(dòng)化研究所則基于先前所設(shè)計(jì)的帶自適應(yīng)先驗(yàn)協(xié)方差的最小二乘法，對每個(gè)網(wǎng)格的邊際對數(shù)似然進(jìn)行了評估，從而確定了2011 年秋季在歐洲探測到的131I 泄漏的源位置和釋放率[23] ?！昂谩钡拇鷥r(jià)函數(shù)的需要滿足：對源參數(shù)敏感，對模型誤差、監(jiān)測誤差和異常值魯棒。一般認(rèn)為平方差對釋放量敏感，而相關(guān)系數(shù)（ corr ）對釋放時(shí)間敏感，因此有學(xué)者通過結(jié)合平方差和相關(guān)系數(shù)，提出了式（4）中的代價(jià)函數(shù)[24] ，成功反演出第三次朝核試驗(yàn)的泄漏源位置和釋放率，與已發(fā)表的結(jié)果相近。

除了代價(jià)函數(shù)的形式，代價(jià)函數(shù)的優(yōu)化方法也非常重要。大部分文獻(xiàn)采用的是基于梯度的方法，比如擬牛頓法[1-2，21-22] ，實(shí)際上也可以采用元啟發(fā)式方法對代價(jià)函數(shù)進(jìn)行優(yōu)化，比如遺傳算法、模擬退火算法和粒子群算法，但是元啟發(fā)式算法不適用于變速釋放行為的泄漏源反演，因?yàn)樵诿總€(gè)時(shí)間步內(nèi)的釋放率都作為源變量，求解空間過大，難以達(dá)到全局最優(yōu)。

2. 2 貝葉斯推斷方法

貝葉斯推斷方法將概率分布引入反演問題中，從而可以解釋輸入數(shù)據(jù)的不確定性。具體地，貝葉斯準(zhǔn)則可以表達(dá)為[18] ：

式中， p（s） 表示源參數(shù)的先驗(yàn)分布，假設(shè)源位置和釋放率的分布是獨(dú)立的，則 p（s） = p（r） p（q） 。

先驗(yàn)分布即寫為這些參數(shù)概率分布的乘積，如果對源參數(shù)信息一無所知，可假設(shè)為均勻分布，由參數(shù)的上下限進(jìn)行約束，源參數(shù)的先驗(yàn)即為常數(shù)。

似然函數(shù) p（μ s） 描述的是源參數(shù) s 的模擬結(jié)果和監(jiān)測的差異，最常見的是高斯似然函數(shù)[22] ：

式中， R 為協(xié)方差矩陣，由模擬誤差和監(jiān)測誤差決定。在計(jì)算時(shí)一般會(huì)對等式（6）取對數(shù)，考慮到監(jiān)測數(shù)據(jù)數(shù)量級往往跨度較大，所以也會(huì)對監(jiān)測數(shù)據(jù)取對數(shù)[25] ，形式如下：

式中， μt 為極小量，可避免零值進(jìn)行對數(shù)計(jì)算時(shí)的失敗，這里是假設(shè)等式（6）中的R = rI ，r 代表觀測-模擬誤差平均協(xié)方差，I 為單位矩陣。

對于變速釋放，如果考慮每個(gè)時(shí)間步的釋放率，參數(shù)維度會(huì)顯著增加，可能導(dǎo)致似然函數(shù)出現(xiàn)多個(gè)局部極值點(diǎn)，即參數(shù)優(yōu)化過程陷入局部最優(yōu)，難以收斂至全局最優(yōu)解。此外，所需的迭代次數(shù)會(huì)呈指數(shù)級增長，這種計(jì)算代價(jià)極大地影響算法的收斂性，尤其在高維空間中更難迅速逼近后驗(yàn)分布的高概率區(qū)域。針對這種情況，通常有兩種處理方法：一種是通過降低時(shí)間分辨率，僅考慮有限時(shí)間步內(nèi)的釋放率[25-26] ;另一種是將變速釋放簡化為矩形釋放[27] 。這些方法可以有效減少計(jì)算復(fù)雜性，提高算法的可行性。

在使用高斯似然且先驗(yàn)為均勻分布的情況下，由于源參數(shù)的先驗(yàn)概率是常數(shù)，所以后驗(yàn)分布與似然函數(shù)同為高斯分布的形式。當(dāng)監(jiān)測數(shù)據(jù)稀疏且?guī)в性肼晻r(shí)，模擬-觀測誤差的方差r 可能變得非常大，這會(huì)使后驗(yàn)分布變得非常平坦和寬廣，接近均勻分布。此外，如果r 的大小未能正確設(shè)定，后驗(yàn)分布將無法準(zhǔn)確反映參數(shù)估計(jì)中的不確定性，從而影響反演結(jié)果的可靠性。因此很多研究都對似然函數(shù)和協(xié)方差參數(shù)的設(shè)計(jì)進(jìn)行了討論。法國國家輻射防護(hù)和核安全研究所利用貝葉斯推論方法對106Ru 泄漏事件進(jìn)行反演，比較了高斯、柯西和拉普拉斯似然函數(shù)對后驗(yàn)分布估計(jì)的結(jié)果，并提出了兩種對監(jiān)測分類并對不同類別分配不同協(xié)方差參數(shù)的方法，而后將各個(gè)協(xié)方差參數(shù)作為未知量估計(jì)后驗(yàn)分布，結(jié)果表明柯西似然函數(shù)的效果最好，且通過對協(xié)方差的建模更準(zhǔn)確地量化了反演結(jié)果的不確定性[25] 。加拿大衛(wèi)生部輻射防護(hù)局則分別構(gòu)建非零監(jiān)測和零值監(jiān)測的似然函數(shù)，并以乘數(shù)的形式表示似然函數(shù)中未考慮到的模型誤差，結(jié)果表明針對不同監(jiān)測考慮不同的模型誤差有助于改善貝葉斯反演的結(jié)果[26] 。先驗(yàn)信息的引入同樣可以改進(jìn)貝葉斯反演的結(jié)果，尤其是對降低求解維度和提高收斂性有很大幫助，但是目前對先驗(yàn)信息的研究較少。

蒙特卡羅（MC）采樣方法用于決定源參數(shù)的后驗(yàn)概率密度分布 p（s μ ） ，基于 MC 方法，大量文獻(xiàn)設(shè)計(jì)了更加高效的采樣方法，比如馬爾可夫鏈蒙特卡羅（MCMC），序列蒙特卡羅方法（SMC）和差分蒙特卡羅（DEMC）[18] 。常見的基于MCMC方法的貝葉斯計(jì)算流程如圖4 所示。

2. 3 源參數(shù)反演方法比較和討論

本質(zhì)上講，迭代優(yōu)化方法是傳統(tǒng)源項(xiàng)反演方法在空間維度上的迭代。迭代優(yōu)化方法只能得到源參數(shù)的單點(diǎn)估計(jì)結(jié)果，而且需要遍歷計(jì)算域內(nèi)所有源位置，因此計(jì)算成本較高。由于對每個(gè)位置都需要求解釋放率，監(jiān)測數(shù)據(jù)過少會(huì)導(dǎo)致釋放率的不確定性變高。但迭代優(yōu)化方法計(jì)算流程簡單，可與傳統(tǒng)源項(xiàng)反演方法無縫銜接，適用于各種釋放行為下的泄漏源參數(shù)反演，可基于此類方法排除一些不切實(shí)際的源參數(shù)組合。

貝葉斯推論方法是優(yōu)化方法的概率化表達(dá)，其計(jì)算成本與源參數(shù)空間大小呈指數(shù)關(guān)系。該方法最大的缺陷是模型的不確定性和觀測的不確定性難以量化，收斂難度大。但貝葉斯推論方法可以利用非常少的監(jiān)測數(shù)據(jù)推斷出源參數(shù)，在最近發(fā)生的一些放射性泄漏事件中得以廣泛應(yīng)用，加拿大衛(wèi)生部輻射防護(hù)局還將該方法設(shè)計(jì)成了應(yīng)用軟件FREAR （ https：/ / gitlab. com/ trDMt2er/FREAR）[27] 。核突發(fā)事件產(chǎn)生的源項(xiàng)釋放特征必然是變速的，當(dāng)氣象條件較為穩(wěn)定時(shí)，變速釋放對預(yù)測濃度的影響主要表現(xiàn)為濃度大小的線性波動(dòng)。在此情況下，即使假設(shè)瞬時(shí)或均勻釋放，仍能有效反映源位置對預(yù)測濃度的影響，因此對位置估計(jì)的準(zhǔn)確性影響較小。位置確定后，可進(jìn)一步采用源項(xiàng)反演方法估計(jì)變速釋放率。然而，當(dāng)氣象條件復(fù)雜多變時(shí)，變速釋放率與時(shí)變氣象場的聯(lián)合效應(yīng)在簡化假設(shè)下被忽略，可能導(dǎo)致源位置參數(shù)的估計(jì)出現(xiàn)不切實(shí)際的結(jié)果。此外，若數(shù)據(jù)極度稀疏，觀測信息無法有效約束解，問題病態(tài)性將變得極其嚴(yán)重，變速釋放假設(shè)會(huì)導(dǎo)致后驗(yàn)分布難以收斂。因此，在實(shí)際使用中，建議基于不同假設(shè)對源參數(shù)進(jìn)行反演，比較各后驗(yàn)分布的收斂性和估計(jì)結(jié)果，從而綜合判斷可能的泄漏源位置。

無論是迭代優(yōu)化方法還是貝葉斯推論方法，往往需要運(yùn)行多次大氣擴(kuò)散模型，計(jì)算代價(jià)非常高，因此有學(xué)者通過擬合機(jī)器學(xué)習(xí)或深度學(xué)習(xí)模型替代擴(kuò)散模型，從而加速反演過程。美國勞倫斯利弗莫爾國家實(shí)驗(yàn)室利用決策樹擬合WRFFLEXPART的輸入變量和輸出指標(biāo)之間的關(guān)系，利用擬合好的模型替代大氣擴(kuò)散模型，極大降低了正向模擬的次數(shù)，并且可以利用決策樹的特征重要性判斷機(jī)制得出輸入變量對輸出結(jié)果的影響排序[28] 。為了提高迭代優(yōu)化方法中源參數(shù)的搜索效率或?yàn)樨惾~斯推論方法提供合理的先驗(yàn)分布，也可以采用快速方法對源位置進(jìn)行初猜，比如利用SRS 場和監(jiān)測之間的相關(guān)性[3] ，或基于反向軌跡的統(tǒng)計(jì)分析[5] 。

3 結(jié)論與建議

放射性泄漏源參數(shù)反演方法在實(shí)際應(yīng)用時(shí)需滿足四個(gè)特性：第一是有效性，需要精確且定量地估計(jì)出源位置和釋放率;第二是快速性，為防止放射性物質(zhì)的擴(kuò)散對公眾和環(huán)境造成巨大威脅，需要快速地估計(jì)出源參數(shù)，及時(shí)進(jìn)行后果評估;第三是可靠性，需要給出反演結(jié)果的不確定性分析;第四是魯棒性，反演方法需要對自身的參數(shù)和監(jiān)測數(shù)據(jù)不敏感。相比之下，貝葉斯推論方法可以利用少量監(jiān)測數(shù)據(jù)提供準(zhǔn)確的源參數(shù)估計(jì)結(jié)果及置信區(qū)間，更滿足實(shí)際應(yīng)用需求，但是貝葉斯推論方法需要選擇合適的似然函數(shù)，且協(xié)方差參數(shù)的選取是至關(guān)重要的，未來需要更多結(jié)合實(shí)際泄漏情景進(jìn)行模擬誤差和監(jiān)測誤差的建模設(shè)計(jì)。除此之外，先驗(yàn)分布的設(shè)計(jì)也可以改進(jìn)貝葉斯推論方法的收斂性，進(jìn)一步提高其準(zhǔn)確性和快速性。