核事故實(shí)時(shí)釋放量集合卡爾曼濾波反演算法研究

唐秀歡，李 華，包利紅

（西北核技術(shù)研究所，陜西 西安 710024）

福島核事故后，根據(jù)核設(shè)施周圍的監(jiān)測(cè)數(shù)據(jù)來(lái)反推［1-2］源項(xiàng)的方法即源項(xiàng)反演越來(lái)越得到重視。目前核事故源項(xiàng)反演技術(shù)有最優(yōu)插值法［3］、遺傳算法［4］、卡爾曼濾波及由其發(fā)展的擴(kuò)展卡爾曼濾波和集合卡爾曼濾波（EnKF）［5-6］等。卡爾曼濾波是一種最優(yōu)遞歸資料處理算法，其中EnKF被認(rèn)為具有較廣的應(yīng)用空間［7］。Astrup等［5-6］采用擴(kuò)展卡爾曼濾波方法跟蹤了單煙團(tuán)的擴(kuò)散過(guò)程，反演了煙團(tuán)位置和煙團(tuán)活度參數(shù)，通過(guò)EnKF 處理，給出了沉降模塊需要的不確定度最佳估計(jì)值。Zheng等［8］結(jié)合隨機(jī)游走擴(kuò)散模型采用EnKF 同化了核素釋放率，實(shí)驗(yàn)中釋放率不確定度減少了80%。文獻(xiàn)［9］探討了卡爾曼濾波反演核事故短時(shí)釋放率的可行性，給出了卡爾曼濾波基本流程和性能，但由于模式誤差的缺失以及線性系統(tǒng)的要求，其應(yīng)用仍受到較大限制。不同的是，EnKF［10］基于隨機(jī)動(dòng)力預(yù)測(cè)理論發(fā)展而來(lái)，用蒙特卡羅短期集合預(yù)報(bào)方法來(lái)估計(jì)預(yù)測(cè)誤差協(xié)方差，具有解決預(yù)測(cè)誤差協(xié)方差矩陣估計(jì)困難及非線性系統(tǒng)近似等問(wèn)題的能力。鑒于此，本文研究適用于核事故現(xiàn)場(chǎng)實(shí)時(shí)釋放量快速反演的EnKF基本算法流程，優(yōu)化EnKF 算法基本參數(shù)，為核事故應(yīng)急源項(xiàng)反演提供新的技術(shù)途徑。

1 非線性源項(xiàng)狀態(tài)模型

對(duì)于已知的核設(shè)施，核事故時(shí)泄漏位置容易觀測(cè)或分析得到，如可能是煙囪、屋頂或安全殼，放射性核素釋放量是核應(yīng)急關(guān)注的重點(diǎn)。核設(shè)施發(fā)生核事故時(shí)，設(shè)核素釋放量狀態(tài)方程如下：

式中：X 為 狀態(tài)向量；k 為 時(shí) 間 段；Mk為 未 知 的變化函數(shù)；qk為過(guò)程激勵(lì)噪聲。

式（1）表明，核素釋放量在各時(shí)間段可能是不同的，以釋放率表示時(shí)，Xk為k 時(shí)間段平均釋放率；以釋放量表示時(shí)，Xk為k 時(shí)間段積分釋放量。

測(cè)量方程為非線性方程：

式中：Zk為源項(xiàng)經(jīng)大氣擴(kuò)散k 時(shí)間段后在核設(shè)施周圍環(huán)境地面空氣的積分活度濃度測(cè)量值；rk為觀測(cè)噪聲；Hk為觀測(cè)矩陣。

式（2）表明，k時(shí)間段前釋放的源項(xiàng)對(duì)k 時(shí)刻的活度濃度依然有影響。

核素大氣擴(kuò)散采用高斯多煙團(tuán)模型計(jì)算［9］，在t時(shí)刻，地面空氣的瞬時(shí)活度濃度與事故源項(xiàng)的釋放率、釋放時(shí)間等有關(guān)，假設(shè)在t=0時(shí)，第i煙團(tuán)釋放出來(lái)（時(shí)段釋放率為Xi），第i煙團(tuán)在tw時(shí)刻即第w 時(shí)段在某網(wǎng)格（x，y，0）上的地面空氣活度濃度為：

2 EnKF的非線性源項(xiàng)跟蹤算法

EnKF反演以核設(shè)施周圍布置的n 個(gè)空間測(cè)量點(diǎn)數(shù)據(jù)為輸入，跟蹤每個(gè)時(shí)間段核素釋放量，流程如下。

1）產(chǎn)生核事故某時(shí)段核素釋放量集合成員數(shù)為N 的預(yù)估計(jì)值集合Xi。

其中：A 為狀態(tài)轉(zhuǎn)移矩陣，假定事故時(shí)釋放量恒定，A=1；下標(biāo)i代表第i個(gè)集合成員，取1，…，N；Xi（n|n-1）為利用上一空間測(cè)量點(diǎn)預(yù)測(cè)的核素釋放量結(jié)果；Xi（n-1|n-1）為上一空間測(cè)量點(diǎn)最優(yōu)的核素釋放量結(jié)果。

2）計(jì)算卡爾曼增益矩陣Kg。EnKF 在計(jì)算增益矩陣時(shí)直接計(jì)算P（n|n-1）·和Hn·P（n|n-1）·。

其中，R=E（rn，rn）為觀測(cè)方程誤差協(xié)方差矩陣。

3）根據(jù)預(yù)測(cè)值和測(cè)量值計(jì)算核事故核素釋放量，更新估計(jì)值X。

4）進(jìn)入下一個(gè)空間測(cè)量點(diǎn)，返回步驟1。

繼續(xù)使用該時(shí)間段上的其他空間測(cè)量點(diǎn)數(shù)據(jù)進(jìn)行反演，最后計(jì)算得到該時(shí)間段上的釋放量最佳估計(jì)值，即為最后一個(gè)空間點(diǎn)的N 個(gè)分析樣本的平均值：

5）進(jìn)入下一時(shí)刻，返回步驟1，連續(xù)跟蹤到最后觀測(cè)值產(chǎn)生的時(shí)刻。

3 反演仿真結(jié)果及分析

3.1 實(shí)驗(yàn)參數(shù)

為了驗(yàn)證EnKF 方法實(shí)時(shí)反演核事故核素釋放量的可行性及其精度，本文設(shè)計(jì)了核事故核素釋放量隨時(shí)間推進(jìn)而變化的反演仿真實(shí)驗(yàn)，釋放量真值變化如下：1）Q=1010+5×109sin（0.314T）；2）Q=1010+109T。數(shù)值實(shí)驗(yàn)設(shè)計(jì)如下：1）假設(shè)釋放高度為35m，釋放時(shí)大氣穩(wěn)定度為D 級(jí)，10m 高度風(fēng)速為4 m／s，風(fēng)向?yàn)镾E，釋放時(shí)間為40min；2）核素監(jiān)測(cè)點(diǎn)放置在核設(shè)施外半徑200、300、400、500、600m處，每個(gè)半徑按π／4 角度平均布置8 個(gè)，共計(jì)40個(gè)，計(jì)算時(shí)假設(shè)只有下風(fēng)向及其鄰近的監(jiān)測(cè)點(diǎn)響應(yīng)核事故的核素釋放，以核設(shè)施為原點(diǎn)，取下風(fēng)向的3π／4扇面監(jiān)測(cè)點(diǎn)數(shù)據(jù)為輸入，設(shè)有效數(shù)據(jù)點(diǎn)為9個(gè)；3）監(jiān)測(cè)數(shù)據(jù)為釋放量真值經(jīng)大氣擴(kuò)散后在地面1 m 處空氣積分活度濃度加入白噪聲后所得，假設(shè)在事故后4 min 時(shí)第一批9個(gè)4 min 積分活度濃度內(nèi)監(jiān)測(cè)數(shù)據(jù)傳回，每隔2min下一批2min積分活度濃度數(shù)據(jù)傳回，每個(gè)數(shù)據(jù)點(diǎn)作一次EnKF 反演，進(jìn)行數(shù)據(jù)處理并反演每一時(shí)段的事故釋放量，最后一個(gè)測(cè)量點(diǎn)數(shù)據(jù)反演處理后作為該時(shí)段的反演結(jié)果；4）R 用過(guò)程噪聲相對(duì)誤差10%來(lái)計(jì)算；5）觀測(cè)時(shí)間段的事故釋放量除以每個(gè)時(shí)間步長(zhǎng)得到時(shí)段釋放率，同時(shí)給出觀測(cè)時(shí)段積分釋放量和釋放率的反演結(jié)果。

3.2 實(shí)驗(yàn)結(jié)果

正弦變化型和直線變化型的源項(xiàng)釋放EnKF反演結(jié)果以及其反演過(guò)程相對(duì)誤差及各時(shí)段最終結(jié)果相對(duì)誤差示于圖1。

由圖1a可見(jiàn)，在第一批觀測(cè)點(diǎn)進(jìn)行EnKF反演后，均可對(duì)積分釋放量和釋放率進(jìn)行跟蹤反演，反演曲線的趨勢(shì)和量級(jí)與真值基本相符。從跟蹤效果上看，積分釋放量的跟蹤總體效果優(yōu)于釋放率的跟蹤，原因在于EnKF 反演使用的模型參數(shù)是積分釋放量。由圖1b可見(jiàn)，在第一批數(shù)據(jù)引入EnKF 反演時(shí)，反演過(guò)程相對(duì)誤差較大，隨后出現(xiàn)周期性變化。圖1c所示的源項(xiàng)隨時(shí)間不斷增大，變化幅度較圖1a的劇烈。從反演結(jié)果看，EnKF 也能對(duì)此類型的釋放進(jìn)行反演，反演曲線的趨勢(shì)和量級(jí)也與真值基本相符。圖1d所示反演過(guò)程的相對(duì)誤差與圖1b所示的趨勢(shì)類似，在整個(gè)反演過(guò)程中均保持同一量級(jí)的鋸齒形；在每個(gè)時(shí)段引入新的隨機(jī)擾動(dòng)后，該時(shí)段反演開(kāi)始時(shí)相對(duì)誤差最大，反演迭代中逐漸變小，符合EnKF 隨機(jī)擾動(dòng)變化的規(guī)律。從圖1b、d看，每個(gè)時(shí)段最終積分誤差較小，基本控制在5%以內(nèi)，有效地收斂到真值。圖1表明，EnKF能對(duì)變化的源項(xiàng)做出反應(yīng)，對(duì)源項(xiàng)目標(biāo)跟蹤具有可行性和有效性。

3.3 集合成員數(shù)、空間測(cè)量點(diǎn)數(shù)對(duì)EnKF的影響

取不同集合成員數(shù)為10、50、100、150、200，以之為變量，結(jié)合不同空間測(cè)量點(diǎn)數(shù)，研究集合成員數(shù)對(duì)EnKF分析的影響，結(jié)果示于圖2。

理論上，EnKF 以基于集合樣本的統(tǒng)計(jì)量來(lái)近似描述真實(shí)的統(tǒng)計(jì)信息，樣本數(shù)的選取對(duì)模型計(jì)算效果的影響是決定性的。根據(jù)大數(shù)定理，隨集合樣本數(shù)的遞增，概率密度函數(shù)的誤差趨于零，收斂階為1／N0.5。由圖2 可見(jiàn)，空間測(cè)量點(diǎn)較多時(shí)（圖2a），在單一時(shí)刻可進(jìn)行多次EnKF反演，其相對(duì)誤差均較小，對(duì)于不同集合數(shù)的反演，所選取的較大樣本數(shù)并未表現(xiàn)出足夠的精度優(yōu)勢(shì)；測(cè)量點(diǎn)較少時(shí)（圖2b），每一時(shí)刻的反演可能不夠充分，相對(duì)誤差變化較為劇烈，集合成員數(shù)為10 時(shí)，反演發(fā)散，而集合數(shù)越大反演效果越好，變化趨勢(shì)符合大數(shù)定理。因此，具有較多的空間測(cè)點(diǎn)時(shí)，取集合數(shù)為50 的EnKF 反演基本取得較佳的反演效果。

圖1 正弦變化型（a，b）和直線變化型（c，d）的源項(xiàng)釋放EnKF反演結(jié)果Fig.1 EnKF inversion results of sine release（a，b）and linear release（c，d）

圖2 不同集合數(shù)對(duì)EnKF反演的影響Fig.2 Influence of different ensemble numbers on EnKF inversion

3.4 初始值對(duì)EnKF的影響

對(duì)初始值的猜測(cè)根據(jù)來(lái)自對(duì)事故本身的監(jiān)測(cè)、安全分析的初步判斷以及測(cè)量點(diǎn)監(jiān)測(cè)數(shù)據(jù)等，為研究初始值對(duì)EnKF 反演的影響，以產(chǎn)生第一批監(jiān)測(cè)數(shù)據(jù)時(shí)刻的釋放量源項(xiàng)模擬真值的1／1 000、1／100、1／10、10、100 倍為初始值，觀察其對(duì)反演結(jié)果的影響，結(jié)果如圖3 所示。由圖3可見(jiàn)，以1／1 000源項(xiàng)真值為初始值時(shí)，EnKF反演不能跟蹤源項(xiàng)的真值，反演所得結(jié)果均遠(yuǎn)小于真值；以1／100源項(xiàng)真值為初始值時(shí)，EnKF反演經(jīng)過(guò)一段時(shí)間后方能跟蹤源項(xiàng)的真值，反演所得結(jié)果偏小；以10、100倍源項(xiàng)真值為初始值時(shí)，EnKF 反演較易追蹤到源項(xiàng)真值。可見(jiàn)，偏大的初始值有利于EnKF 反演。在實(shí)際應(yīng)用中，應(yīng)根據(jù)現(xiàn)場(chǎng)的其他觀測(cè)結(jié)果以及對(duì)事故的認(rèn)識(shí)，合理估計(jì)偏大的釋放初始源項(xiàng)，為EnKF反演提供恰當(dāng)?shù)某跏贾怠?/p>

圖3 不同初始值對(duì)EnKF反演的影響Fig.3 Influence of different initial values on EnKF inversion

3.5 擾動(dòng)值、觀測(cè)誤差對(duì)EnKF的影響

固定測(cè)量相對(duì)誤差為10%，采用添加隨機(jī)擾動(dòng)法對(duì)每一時(shí)段初期的集合成員進(jìn)行一次擾動(dòng)，擾動(dòng)值取上時(shí)刻均值的10、5、2、1、0.5、0.1及0倍，觀察擾動(dòng)值對(duì)反演結(jié)果的影響。固定擾動(dòng)值為10倍上時(shí)刻均值，取測(cè)量相對(duì)誤差的1 000%、100%、50%、20%、10%、1%及0.1%作為觀測(cè)誤差，觀察觀測(cè)誤差對(duì)反演結(jié)果的影響，結(jié)果如圖4 所示。由圖4a 可見(jiàn)，無(wú)擾動(dòng)（0倍）、較小擾動(dòng)如0.1倍、0.5倍均值，反演相對(duì)誤差逐漸增大，反演值趨向于發(fā)散；1倍均值擾動(dòng)反演相對(duì)誤差保持在10%左右，2、5、10倍擾動(dòng)反演相對(duì)誤差總體上基本控制在10%內(nèi)，但2 倍均值擾動(dòng)造成的反演相對(duì)誤差總體較5倍和10倍的大，5倍均值反演過(guò)程相對(duì)誤差局部波動(dòng)較10 倍的大，10 倍均值擾動(dòng)反演相對(duì)誤差最小，即較大擾動(dòng)值的反演效果較佳。在觀測(cè)相對(duì)誤差影響方面，由圖4b可見(jiàn)，100%以上相對(duì)誤差造成反演迭代的發(fā)散，20%～50%相對(duì)誤差造成的反演誤差約在10%以內(nèi)，10%相對(duì)誤差造成的反演誤差則在5%以內(nèi)，觀測(cè)相對(duì)誤差越小則反演精度越高。

3.6 討論

源項(xiàng)反演還有一些其他的優(yōu)化技術(shù)，如遺傳算法、網(wǎng)絡(luò)神經(jīng)模擬等［7］，反演的參數(shù)、原理、效率也不盡相同，EnKF 的優(yōu)勢(shì)在于可實(shí)時(shí)跟蹤源項(xiàng)釋放量的變化，較快地完成反演過(guò)程。核事故源項(xiàng)釋放參數(shù)需要反演的可能有釋放量、釋放高度、釋放時(shí)間等，就核事故形成的巨大危害輻射場(chǎng)而言，釋放量及其隨時(shí)間的變化是問(wèn)題的核心。對(duì)于其他反演量，如釋放高度、釋放時(shí)間，由于它們之間的相互關(guān)系未知，因此多個(gè)參數(shù)的同時(shí)反演存在技術(shù)困難，目前只有單個(gè)參數(shù)的反演是可行的。

圖4 不同擾動(dòng)值（a）及不同觀測(cè)相對(duì)誤差（b）對(duì)EnKF反演的影響Fig.4 Influence of different perturbations（a）and observation relative errors（b）on EnKF

離散卡爾曼濾波適用于線性的源項(xiàng)跟蹤模型，只適合短期釋放的平均釋放類型，其反演需要測(cè)量誤差和模型誤差作為輸入?yún)?shù)［9］。EnKF可用于非線性源項(xiàng)釋放系統(tǒng)，可根據(jù)源項(xiàng)釋放變化情況追蹤其變化率，反演結(jié)果雖受模型誤差的影響，但無(wú)需該參數(shù)作為輸入，比卡爾曼線性濾波具有優(yōu)越性。本文反演仿真中釋放源項(xiàng)是虛擬設(shè)計(jì)的，大氣高斯多煙團(tuán)擴(kuò)散模型是基于正態(tài)分布簡(jiǎn)化的，而在實(shí)際應(yīng)用中會(huì)存在一系列的計(jì)算誤差，如風(fēng)場(chǎng)不確定性、沉降和化學(xué)反應(yīng)過(guò)程近似及模式分辨率等。另外，隨時(shí)間的延長(zhǎng)，源項(xiàng)擴(kuò)散的煙團(tuán)數(shù)目也逐漸增多，導(dǎo)致了反演計(jì)算量的增大。為保證反演盡可能覆蓋長(zhǎng)的時(shí)間段，建議對(duì)煙團(tuán)時(shí)間長(zhǎng)度進(jìn)行放大，或?qū)σ恍┻h(yuǎn)離固定位置的煙團(tuán)進(jìn)行記錄和過(guò)濾，只計(jì)算對(duì)近期釋放的活度濃度貢獻(xiàn)較大的煙團(tuán)。目前，在核事故現(xiàn)場(chǎng)可實(shí)時(shí)測(cè)量的數(shù)據(jù)多為地面空氣γ吸收劑量率，由于核事故釋放源項(xiàng)核素多樣、化學(xué)成分復(fù)雜，每種核事故具有不同的核素釋放譜，核素發(fā)射的射線受現(xiàn)場(chǎng)空間分布的影響，散射、沉積于空氣中的能量多樣化，地面空氣γ吸收劑量率與核事故現(xiàn)場(chǎng)復(fù)雜核素的活度濃度還存在著復(fù)雜未知的關(guān)系。本文采用地面空氣積分活度濃度作為測(cè)量數(shù)據(jù)輸入，只是在原理上證實(shí)了EnKF 技術(shù)的可行性和有效性，現(xiàn)場(chǎng)地面空氣γ吸收劑量率測(cè)量值的輸入問(wèn)題可作為未來(lái)研究的重點(diǎn)。

4 結(jié)論

本文在源項(xiàng)反演跟蹤系統(tǒng)中引入了EnKF技術(shù)，建立了與高斯多煙團(tuán)擴(kuò)散模型相結(jié)合的核事故實(shí)時(shí)釋放量EnKF 基本算法流程。結(jié)果顯示，集合成員數(shù)為50、偏大初始值以及較大擾動(dòng)值為輸入的反演效果較為理想，優(yōu)化后的流程在虛擬數(shù)據(jù)反演中積分誤差在5%以內(nèi)。研究表明，核事故實(shí)時(shí)釋放量EnKF 算法具有較好的有效性和可行性，原則上可應(yīng)用于核事故源項(xiàng)釋放量連續(xù)實(shí)時(shí)跟蹤反演。

［1］ 陳曉秋，楊端節(jié)，李冰.利用福島第一核電廠事故期間環(huán)境監(jiān)測(cè)資料反推事故釋放源項(xiàng)［J］.核化學(xué)與放射化學(xué)，2012，34（2）：83-87.CHEN Xiaoqiu，YANG Duanjie，LI Bing.Reverse estimation of accidental release amounts from Fukushima Daiichi Nuclear Power Plant by environmental monitoring data［J］.Journal of Nuclear and Radiochemistry，2012，34（2）：83-87（in Chinese）.

［2］ 程衛(wèi)亞，楊宏偉，陳凌，等.利用航測(cè)數(shù)據(jù)反推福島核事故137Cs的釋放量［J］.原子能科學(xué)技術(shù)，2012，46（2）：252-256.CHENG Weiya，YANG Hongwei，CHEN Ling，et al.Deducing total released activity of137Cs in Fukushima’s nuclear accident by results of airborne monitoring［J］.Atomic Energy Science and Technology，2012，46（2）：252-256（in Chinese）.

［3］ JEONG H J，KIM E H，SUH K S，et al.Determination of the source rate released into the environment from a nuclear power plant［J］.Radiation Protection Dosimetry，2005，113（3）：308-313.

［4］ 寧莎莎，蒯琳萍.混合遺傳算法在核事故源項(xiàng)反演中的應(yīng)用［J］.原子能科學(xué)技術(shù)，2012，46（增刊）：469-472.NING Shasha，KUAI Linping.Back calculation of source terms by hybrid genetic algorithm in nuclear power plant accident［J］.Atomic Energy Science and Technology，2012，46（Suppl.）：469-472（in Chinese）.

［5］ ASTRUP P，TURCANU C，PUCH R O，et al.Data assimilation in the early phase：Kalman filtering RIMPUFF，RODOS（RA5）-TN（04）-01［R］.Denmark：Ris?National Laboratory Information Service Department，2004.

［6］ PUCH R O，ASTRUP P.An extended Kalman filter methodology for the plume phase of a nuclear accident，RODOS（RA5）-TN（01）-07［R］.Denmark：Ris?National Laboratory Information Service Department，2002.

［7］ 徐志新，奚樹(shù)人，曲靜原.核事故源項(xiàng)反演技術(shù)及其研究現(xiàn)狀［J］.科技導(dǎo)報(bào)，2007，25（5）：16-20.XU Zhixin，XI Shuren，QU Jingyuan.Review on source inversion technology in analyzing nuclear accidents［J］.Science ＆ Technology Review，2007，25（5）：16-20（in Chinese）.

［8］ ZHENG D Q，LEUNG J K C，LEE B Y，et al.Data assimilation in the atmospheric dispersion model for nuclear accident assessments［J］.Atmospheric Environment，2007，41：2 438-2 446.

［9］ 唐秀歡，包利紅，李華，等.卡爾曼線性濾波反演核設(shè)施核事故中核素釋放率的研究［J］.原子能科學(xué)技術(shù)，2014，48（10）：1 847-1 852.TANG Xiuhuan，BAO Lihong，LI Hua，et al.Radionuclide release rate inversion of nuclear accidents in nuclear facility based on Kalman filter［J］.Atomic Energy Science and Technology，2014，48（10）：1 847-1 852（in Chinese）.

［10］EVENSEN G.The ensemble Kalman filter：Theoretical formulation and practical implementation［J］.Ocean Dynamics，2003，53：343-367.