唐秀歡，李 華，包利紅

（西北核技術研究所，陜西 西安 710024）

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

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

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

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

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

測量方程為非線性方程：

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

式（2）表明，k時間段前釋放的源項對k 時刻的活度濃度依然有影響。

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

2 EnKF的非線性源項跟蹤算法

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

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

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

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

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

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

4）進入下一個空間測量點，返回步驟1。

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

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

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

3.1 實驗參數(shù)

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

3.2 實驗結(jié)果

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

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

3.3 集合成員數(shù)、空間測量點數(shù)對EnKF的影響

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

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

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

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

3.4 初始值對EnKF的影響

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

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

3.5 擾動值、觀測誤差對EnKF的影響

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

3.6 討論

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

圖4 不同擾動值（a）及不同觀測相對誤差（b）對EnKF反演的影響Fig.4 Influence of different perturbations（a）and observation relative errors（b）on EnKF

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

4 結(jié)論

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

［1］ 陳曉秋，楊端節(jié)，李冰.利用福島第一核電廠事故期間環(huán)境監(jiān)測資料反推事故釋放源項［J］.核化學與放射化學，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)亞，楊宏偉，陳凌，等.利用航測數(shù)據(jù)反推福島核事故137Cs的釋放量［J］.原子能科學技術，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］ 寧莎莎，蒯琳萍.混合遺傳算法在核事故源項反演中的應用［J］.原子能科學技術，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］ 徐志新，奚樹人，曲靜原.核事故源項反演技術及其研究現(xiàn)狀［J］.科技導報，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］ 唐秀歡，包利紅，李華，等.卡爾曼線性濾波反演核設施核事故中核素釋放率的研究［J］.原子能科學技術，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.