唐秀歡，包利紅，李 華，萬俊生

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

估算核設(shè)施核事故源項一般有3種方法：監(jiān)測手段、安全分析中的事故分析及事故中核設(shè)施工況的詳細(xì)分析[1]。后兩種方法往往過于保守，或在事故時難以實(shí)施，例如福島核事故早期全廠斷電，詳細(xì)分析難以完成。因此，根據(jù)核設(shè)施周圍的監(jiān)測數(shù)據(jù)來反推源項的方法即源項反演得到重視。目前核事故源項反演主要分為以下4類[2]：最優(yōu)插值法、遺傳算法、人工神經(jīng)網(wǎng)絡(luò)、卡爾曼濾波及由其發(fā)展的擴(kuò)展卡爾曼濾波和集合卡爾曼濾波。Jeong等[3]使用最優(yōu)插值法估計了Yeoung-Kwang核設(shè)施示蹤氣體實(shí)驗的點(diǎn)源釋放率。陳曉秋等[4]采用線性化處理反推了福島核事故的核素釋放量。寧沙沙等[5]通過遺傳編碼搜索了最佳的核素釋放率、釋放高度及釋放位置參數(shù)；陳媛華[6]建立了環(huán)境污染事件源項反演的遺傳算法程序；陳軍明等[7]使用遺傳算法和多點(diǎn)源擴(kuò)散模式，從控制點(diǎn)濃度成功地反演了點(diǎn)源排放清單。Astrup等[8-9]采用擴(kuò)展卡爾曼方法跟蹤了單煙團(tuán)的擴(kuò)散過程，反演了煙團(tuán)位置和煙團(tuán)活度參數(shù)。朱江等[10]、Zheng等[11]研究了集合卡爾曼方法對污染源的反演應(yīng)用，討論了污染源反演問題的統(tǒng)計方法。

核事故應(yīng)急監(jiān)測數(shù)據(jù)往往是核設(shè)施周邊環(huán)境固定、移動的劑量或地面空氣活度濃度，以監(jiān)測數(shù)據(jù)反演源項，必然涉及測量和預(yù)測的誤差；同時，由于監(jiān)測數(shù)據(jù)具有實(shí)時和連續(xù)性，反演過程應(yīng)具有連續(xù)數(shù)據(jù)遞歸處理的能力?？柭鼮V波是一種最優(yōu)遞歸資料處理算法，它綜合利用一切可能的觀測信息以及模式和觀測信息的誤差統(tǒng)計特征，使估計的統(tǒng)計誤差達(dá)到最小。本文研究卡爾曼濾波反演核事故核素釋放率算法，實(shí)現(xiàn)其在核設(shè)施核事故應(yīng)急中的有效應(yīng)用。

1 原理和方法

1.1 狀態(tài)和測量方程

對于非移動的核設(shè)施，核應(yīng)急重點(diǎn)關(guān)注的核素釋放率狀態(tài)方程如下：

Xk+1=Xk+qk

(1)

式中：X為狀態(tài)向量；k為采樣系列；qk為過程激勵噪聲。

測量方程為線性方程：

Zk=HkXk+rk

(2)

式中：Zk為源項經(jīng)大氣擴(kuò)散后在核設(shè)施周圍環(huán)境地面空氣的瞬時活度濃度測量值；rk為觀測噪聲；Hk為觀測矩陣。采用文獻(xiàn)[12]描述的高斯多煙團(tuán)模型計算放射性核素的大氣擴(kuò)散，在t時刻，地面空氣的瞬時活度濃度與事故源項的釋放率、釋放時間等有關(guān)，假設(shè)在t=0時，第i煙團(tuán)釋放出來(時段釋放率為Xi)，第i煙團(tuán)在tw時刻即第w時段在某網(wǎng)格(x，y，0)上的地面空氣活度濃度為：

(3)

(4)

(5)

擴(kuò)散參數(shù)計算采用Briggs修正公式[12]，考慮干沉積過程造成的煙羽耗損和大氣混合層對煙羽擴(kuò)散的影響。

1.2 算法流程實(shí)現(xiàn)

將核事故核素釋放率狀態(tài)噪聲和地面活度濃度的觀測噪聲均視為高斯白噪聲，狀態(tài)方程和觀測方程的誤差協(xié)方差矩陣分別為Q和R，則卡爾曼濾波算法步驟如下。

1) 計算核事故核素釋放率預(yù)估計值X：

X(k|k-1)=A·X(k-1|k-1)

(6)

其中：A為狀態(tài)轉(zhuǎn)移矩陣，假定釋放時釋放率恒定，A=1；X(k|k-1)為利用上一狀態(tài)預(yù)測的核素釋放率結(jié)果；X(k-1|k-1)為上一狀態(tài)最優(yōu)的核素釋放率結(jié)果。

2) 計算預(yù)估計協(xié)方差矩陣P：

P(k|k-1)=A·P(k-1|k-1)·A′+Q

(7)

其中，Q=E(qk，qk)為狀態(tài)方程誤差協(xié)方差矩陣，不同時刻的觀測結(jié)果其平行性和準(zhǔn)確性均不同，在迭代中其值亦不同。

3) 計算卡爾曼增益矩陣Kg：

(8)

其中，R=E(rk，rk)為觀測方程誤差協(xié)方差矩陣，當(dāng)測量噪聲極大時，殘余中可用的信息量極低，增益很小，卡爾曼濾波器無法發(fā)揮作用，此時應(yīng)使用先驗估計值。在大氣擴(kuò)散計算中Hk是變量，每次測量數(shù)值的激勵系數(shù)均不同。

4) 根據(jù)預(yù)測值和測量值計算核事故核素釋放率更新估計值X：

X(k|k)=X(k|k-1)+Kg(k)·

(Z(k)-Hk·X(k|k-1))

(9)

5) 計算更新后估計值的協(xié)方差矩陣P：

P(k|k)=(I-Kg(k)·Hk)·P(k|k-1)

(10)

綜合以上各式，可得卡爾曼濾波的計算流程如圖1所示。

先驗估計值采用源項經(jīng)驗判斷結(jié)果，或用最優(yōu)插值法計算第1批次數(shù)據(jù)作為先驗估計值。

圖1 卡爾曼濾波算法流程

2 數(shù)值實(shí)驗與結(jié)果分析

2.1 數(shù)值實(shí)驗

為驗證卡爾曼濾波方法實(shí)時反演核事故核素釋放率的可行性及其精度，本文設(shè)計了以下兩種方案的反演數(shù)值實(shí)驗。第1種，假定核事故核素釋放率是穩(wěn)定均勻的，即不考慮核事故核素釋放率隨時間的變化；第2種，考慮核事故核素釋放率隨采樣點(diǎn)數(shù)而變化。實(shí)驗過程設(shè)計如下：1) 假設(shè)釋放高度為35 m，釋放時大氣穩(wěn)定度為D級，10 m高度風(fēng)速為4 m/s，風(fēng)向為SE，釋放時間為10 min；2) 核素監(jiān)測點(diǎn)放置在核設(shè)施外半徑200、300、400、500、600 m處，每個半徑按π/4角度平均布置8個，共計40個，計算時假設(shè)只有下風(fēng)向及其鄰近的監(jiān)測點(diǎn)響應(yīng)核事故的核素釋放，以核設(shè)施為原點(diǎn)，取包絡(luò)下風(fēng)向的π扇面監(jiān)測點(diǎn)數(shù)據(jù)為輸入，共計25個；3) 以真值加入白噪聲后所得數(shù)據(jù)為監(jiān)測數(shù)據(jù)，在事故后2 min時第1批25個監(jiān)測數(shù)據(jù)傳回，每隔2 min下一批數(shù)據(jù)傳回，每個數(shù)據(jù)點(diǎn)作一次反演，由卡爾曼濾波器進(jìn)行數(shù)據(jù)處理并反演事故釋放率；4) 狀態(tài)方程和觀測方程的誤差協(xié)方差矩陣Q和R用過程噪聲相對誤差10%和觀測相對誤差50%來計算。

2.2 釋放率不隨時間變化條件

在核設(shè)施發(fā)生核事故，核素釋放率不隨時間變化，其為一常數(shù)時，設(shè)核素釋放率為105Bq/s。采用環(huán)境監(jiān)測數(shù)據(jù)及不同初始X0值對核素釋放率參數(shù)進(jìn)行卡爾曼濾波反演，2～10 min內(nèi)真值、測量值及卡爾曼濾波反演結(jié)果示于圖2。

由圖2可見，核設(shè)施核事故核素釋放率不隨時間變化時，卡爾曼濾波法可跟蹤核素釋放率，其反演結(jié)果較測量數(shù)據(jù)平滑。在卡爾曼濾波迭代初始值較真值低時，卡爾曼濾波法反演約10次后追蹤到真值，隨反演點(diǎn)數(shù)的加大，反演標(biāo)準(zhǔn)差逐漸比測量標(biāo)準(zhǔn)差小(圖2a、b)；在卡爾曼濾波迭代初始值較真值高(真值的105倍)時，卡爾曼濾波法經(jīng)60次反演后追蹤到真值，由于初始誤差較大，總體反演標(biāo)準(zhǔn)差比測量標(biāo)準(zhǔn)差高(圖2c、d)；由反演效果看，較小的初始X0值具有較快的追蹤速度和較小的標(biāo)準(zhǔn)差(圖2a、b)。可見，當(dāng)核設(shè)施核事故核素釋放率不隨時間變化時，根據(jù)環(huán)境監(jiān)測數(shù)據(jù)進(jìn)行的卡爾曼濾波法很快收斂，可有效跟蹤核素釋放率。

2.3 釋放率隨采樣點(diǎn)數(shù)變化條件

核設(shè)施發(fā)生核事故時，核素釋放率往往隨時間發(fā)生變化，先假設(shè)核素釋放率真值隨反演點(diǎn)數(shù)線性變化，即X=105t(t=1，2，…，125)；再假設(shè)核素釋放率真值隨點(diǎn)數(shù)非線性變化，即X=|105-102(t-10)2+103t|(t=1，2，…，125)。采用環(huán)境監(jiān)測數(shù)據(jù)及不同初始P0值對核素釋放率參數(shù)進(jìn)行卡爾曼濾波反演，真值、測量值及卡爾曼濾波反演結(jié)果示于圖3。

由圖3可見：在核素釋放率真值線性變化條件下，根據(jù)環(huán)境監(jiān)測數(shù)據(jù)進(jìn)行卡爾曼濾波反演可跟蹤真值的線性變化；同樣在核素釋放率真值非線性變化條件下，根據(jù)環(huán)境監(jiān)測數(shù)據(jù)進(jìn)行卡爾曼濾波反演也能追蹤核素釋放率的變化。在線性和非線性兩種變化條件下，卡爾曼濾波均經(jīng)過約10次反演后即可追蹤到核素釋放率真值。在誤差方面，不同初始P0值對反演結(jié)果影響很小，經(jīng)過約20次反演后卡爾曼濾波的標(biāo)準(zhǔn)差比測量值的標(biāo)準(zhǔn)差小，總體標(biāo)準(zhǔn)差隨真值的增大而增大?？梢姡?dāng)核設(shè)施核事故核素釋放率隨時間變化時，根據(jù)環(huán)境監(jiān)測數(shù)據(jù)進(jìn)行的卡爾曼濾波法約10次追蹤到真值，可有效跟蹤核素釋放率的線性或非線性變化。

2.4 討論

核設(shè)施核事故放射性核素排放可能出現(xiàn)很復(fù)雜的情況：時間上不連續(xù)性，如短時連續(xù)釋放、長時間間斷釋放；空間上不均勻性，如有效釋放高度變化、氣象變化造成的紊流擴(kuò)散等。核事故源項的狀態(tài)向量可能需要包含釋放位置、釋放高度、釋放率等。本文設(shè)計的卡爾曼濾波法主要對線性系統(tǒng)進(jìn)行分析，適用于固定高度和位置、短時連續(xù)排放的核事故核素釋放率參數(shù)反演，且前提條件是假設(shè)環(huán)境監(jiān)測傳回的數(shù)據(jù)經(jīng)過預(yù)先處理后形成核素釋放率數(shù)據(jù)的測量值。常用源項反演最佳插值法是卡爾曼濾波法的簡化版，主要用于線性系統(tǒng)的單變量分析，但對連續(xù)批量數(shù)據(jù)的處理能力不足。實(shí)際情況是釋放高度、釋放位置等多參數(shù)向量的觀測方程和狀態(tài)轉(zhuǎn)移方程均具有很強(qiáng)的非線性。對于非線性特性，可采用擴(kuò)展卡爾曼濾波進(jìn)行基于泰勒級數(shù)展開的線性化處理，得到包括一階近似項的狀態(tài)方程和測量方程的近似表述形式。本文在不同釋放高度位置對釋放高度參數(shù)展開泰勒級數(shù)，分析一階近似項截斷誤差，結(jié)果示于圖4。由圖4可見，在不同展開高度，當(dāng)釋放高度真值不同時，截斷相對誤差不同，真值距展開點(diǎn)越遠(yuǎn)，截斷相對誤差越大，且數(shù)量級可達(dá)1015，遠(yuǎn)超出了反演精度要求。

圖2 釋放率不隨時間變化時卡爾曼濾波反演結(jié)果

圖3 釋放率隨時間變化時卡爾曼濾波反演結(jié)果

圖4 釋放高度擴(kuò)展卡爾曼濾波截斷相對誤差

大多核事故的位置地點(diǎn)是已知的，本文在狀態(tài)向量選取中只考慮釋放率，文獻(xiàn)[8-9]的狀態(tài)向量則更加關(guān)注移動的核事故釋放源，考慮的源項參數(shù)有位置、釋放率、高度等，例如單煙團(tuán)反演系統(tǒng)狀態(tài)向量(x，y，z，q)[8]。對于源項的其他參數(shù)如泄漏時間、泄漏總量等特征量尚未見反演的研究報道。另外，本文測量方程中測量值為核素地面空氣的瞬時活度濃度，而非容易獲得的地面空氣的γ吸收劑量率。原因在于γ空氣吸收劑量率隨核素組成、濃度大小、煙云空間分布等變化而變化，如何通過γ空氣吸收劑量率測量數(shù)據(jù)來確定每個核素的含量仍是復(fù)雜的技術(shù)問題。對于卡爾曼濾波法，未來的研究重點(diǎn)是協(xié)方差數(shù)據(jù)和γ空氣吸收劑量率與測量方程的關(guān)系，特別是核事故應(yīng)急測量的活度以及大氣擴(kuò)散模型的誤差，關(guān)鍵核素的γ空氣吸收劑量率與觀測方程的物理關(guān)系。通過減少模型誤差，簡化輸入?yún)?shù)，從而提高卡爾曼濾波法核素釋放率的反演能力。

3 結(jié)論

根據(jù)核設(shè)施核事故環(huán)境監(jiān)測數(shù)據(jù)的實(shí)時性和連續(xù)性特點(diǎn)，采用與高斯多煙團(tuán)擴(kuò)散模型結(jié)合的卡爾曼濾波法可用于固定高度和位置、短時連續(xù)排放的核事故核素釋放率參數(shù)的反演，這為核設(shè)施核事故應(yīng)急提供一種新的源項反演手段。由于釋放高度參數(shù)的源項反演觀測方程和狀態(tài)轉(zhuǎn)移方程具有較強(qiáng)的非線性，近似處理后截斷誤差不能滿足反演要求，基于近似處理的擴(kuò)展卡爾曼濾波法在釋放高度參數(shù)源項反演的應(yīng)用中有著一定的局限性。本文重點(diǎn)關(guān)注核事故核素釋放率參數(shù)，在將來的工作中，卡爾曼濾波法觀測方程和狀態(tài)轉(zhuǎn)移方程實(shí)際協(xié)方差數(shù)據(jù)的獲取和分析，以及γ空氣吸收劑量率與測量方程物理關(guān)系的研究將是決定該方法能否得到實(shí)際應(yīng)用的關(guān)鍵。

參考文獻(xiàn)：

[1] 王醒宇，康凌. 核事故后果評價方法及其新發(fā)展[M]. 北京：原子能出版社，2003：11.

[2] 徐志新，奚樹人，曲靜原. 核事故源項反演技術(shù)及其研究現(xiàn)狀[J]. 科技導(dǎ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).

[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] 陳曉秋，楊端節(jié)，李冰. 利用福島第一核電廠事故期間環(huán)境監(jiān)測資料反推事故釋放源項[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 environment monitoring date[J]. Journal of Nuclear and Radiochemistry, 2012, 34(2): 83-87(in Chinese).

[5] 寧莎莎，蒯琳萍. 混合遺傳算法在核事故源項反演中的應(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).

[6] 陳媛華. 河流突發(fā)環(huán)境污染事件源項反演及程序設(shè)計[D]. 哈爾濱：哈爾濱工業(yè)大學(xué)，2011.

[7] 陳軍明，徐大海，朱蓉. 遺傳算法在點(diǎn)源擴(kuò)散濃度反演排放源強(qiáng)中的應(yīng)用[J]. 氣象，2002，28(9)：12-16.

CHEN Junming, XU Dahai, ZHU Rong. Application of genetic algorithms to point-source inversion[J]. Meteorological Monthly, 2002, 28(9): 12-16(in Chinese).

[8] 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.

[9] 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.

[10] 朱江，汪萍. 集合卡爾曼平滑和集合卡爾曼濾波在污染源反演中的應(yīng)用[J]. 大氣科學(xué)，2006，30(5)：871-882.

ZHU Jiang, WANG Ping. Ensemble Kalman smoother and ensemble Kalman filter approaches to the joint air quality state and emission estimation problem[J]. Chinese Journal of Atmospheric Sciences, 2006, 30(5): 871-882(in Chinese).

[11] 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.

[12] 胡二邦，陳家宜. 核電廠大氣擴(kuò)散及其環(huán)境影響評價[M]. 北京：原子能出版社，1999.